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SECTION  I 


INTRODUCTION 

The  increasing  use  of'  composite  materials  in  structural  applications,  such  as 
automobiles,  aircraft  and  space  structures,  is  characterized  by  their  high  strength 
(stinness)-tO'Weight  ratio,  low  maintenance  costs  and  the  flexibility  in  tailoring  the 
stiffness  and  strength  to  design  requirements.  As  fiber  reinforced  laminates  h.ive 
played  a  more  important  role  in  high  performance  structures  for  the  last  2  decades,  the 
need  to  have  accurate  stress  and  failure  analysis  become  apparent  for  design  or  repair 
purpose. 

Recent  development  in  the  analysis  of  composite  laminate  coupons  under  uniform 
extension  indicated  that  the  high  interlaminar  stresses  near  the  free-edge  are  mainly 
responsible  for  delamination  failure  [l].  Before  delamination  can  be  predicted  on  the 
basis  of  a  stress-based  failure  criterion,  it  is  essential  that  a  highly  reliable  estimate  of 
interlaminar  stresses  be  available  for  the  given  situation.  However,  it  has  been 
difficult  to  obtain  solutions  for  the  stress  field  because  of  the  anisotropy  as  well  as 
heterogeneity  of  the  material,  and  the  difficulty  of  satisfying  traction-free  boundary 
condition  in  a  solution  procedure  based  on  the  displacement  formulation. 

Consideraole  re.search  efforts  have  been  devoted  to  the  study  of  such  free-edge 
delamination  problem.  These  can  be  classified  as  analytical  and  numerical  approaches. 
The  analytical  solutions  are  based  upon  simple  elastic  approximation  [2,3],  modified 
higher  order  theory  [4],  Galerkin  method  [5],  Perturbation  technique  [6],  Boundary  layer 


theory  [7],  Reissner's  variational  principle  [8,9],  Cjlohal-](x;a]  model  flO]  etc.,  while  the 
numerical  solutions  are  based  on  finite  difference  [ll,12]  and  finite  element  methods 
including  displacement  [13-16],  stress  [l7]  and  hybrid  [18,19]  formulations.  It  was 
found  that  some  of  the  solution  techniques  were  only  applicable  under  certain 
conditions.  For  this  reason,  a  complete  stress  distribution  was  usually  hard  to  obtain. 
Although  results  calculated  from  various  approaches  have  demonstrated  similarities  in 
some  cases,  discrepancies  dn  exist  in  the  magnitude  as  well  as  sign  of  the  computed 
interlaminar  stresses  near  the  free-edge  of  laminate  coupons.  One  example  is  shown  in 
Figure  (1)  in  which  significant  difference  was  obserxed  for  a,  stress  distribution  along 
the  interface  of  [45/-451.  laminate  based  on  various  solution  techniques  [19]. 
Apparently,  one  possible  source  of  these  discrepancies  is  that,  in  these  methods,  the 
continuity  conditions  for  displacements  and  tractions  across  laminate  interfaces  along 
with  traction-free  boundary  condition  along  free-edges  characteristic  of  the  real  life 
situation,  can  only  be  approximated  to  a  limited  extent.  However,  the  credibility  of 
various  methods  in  predicting  the  cr^  distribution  shown  in  Figure  (l)  will  be  judged 
later. 

Due  to  the  presence  of  singular  interlaminar  stresses  near  the  laminate 
free-boundary,  edge  delamination  ass(x;iaied  with  various  types  of  damage  modes,  such 
as  fiber  breakage,  matrix  cracking,  fiber-matrix  debonding,  etc.,  are  observed  to  occur 
under  incremental  loading.  Delamination  can  be  simply  interpreted  as  separation  of 
laminae  from  each  other  in  the  laminate,  and  can  occur  under  static,  impact  or  fatigue 
loading  conditions.  For  the  case  of  a  symmetric  laminate  under  inplane  loading,  the 
strain  components  are  essentially  uniform  throughout  the  laminate.  Due  to  the 
free-edge  effect  the  out-of-plane  interlaminar  stresses,  however,  may  be  sufficiently 
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large  to  damage  the  matrix  material,  which  bonds  adjacent  plies  together,  and  cause 
delamination. 

Two  generic  approaches  are  available  for  investigating  damage  modes  in  composite 
materials.  The  first  approach  involves  a  detailed  stress  analysis  used  in  conjunction 
with  a  failure  criterion  to  predict,  and  measure  experimentally,  the  onset  of  fiber 
fracture,  matrix  cracking  and  delamination.  This  can  be  referred  to  as  the  strength 
characterization  approacli.  In  the  second  approach,  classical  linear  elastic  fracture 

mechanics  can  be  applied  to  characterize  matrix  cracking  and  the  delamination  process. 
Delamination  has  usually  been  isolated  from  the  other  damage  modes  and  treated  as  a 
stable  crack  growth  [20-25],  and  the  basic  character  of  the  strain  energy  release  rate 
has  been  widely  used  to  predict  the  kinematic  behavior  of  delamination.  However, 
experiments  [26]  have  indicated  that  delamination  usually  does  not  produce  a  clean 
surface  between  the  adjacent  plies;  instead  it  is  associated  with  other  types  of  damage 
such  as  matrix  cracking  and  fiber  breakage.  Thus,  use  of  linear  elastic  fracture 
mechanics  approach  to  study  delamination  growth  seems  to  be  inappropriate. 
Meanwhile,  due  to  the  irregular  occurrence  of  various  damage  modes  in  the  form  of 

different  cracking  patterns,  use  of  anisotropic  strength  and  failure  criteria  is  apparently 
superior  to  the  fracture  mechanic  approach  for  the  determination  of  damage 

characteristics  such  as  type  of  failure  mode,  damage  zone  and  crack  growth  behavior 

including  delamination. 

The  primary  objective  of  the  present  research  was  to  develop  a  finite  element 
model  with  a  sound  theoretical  background,  which  could  accurately  and  efficiently 
predict  the  complete  stress  field  of  the  free-edge  stress  problem  in  composite  laminates 
without  resorting  to  any  special  singularity  elements.  The  next  was  to  incorporate 
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various  commonly  known  macroscopic  failure  criteria  into  the  finite  element 
computational  procedure  to  evaluate  the  performance  of  various  criteria  on  the 
determination  of  onset  of  matrix  cracking  and  delamination  in  the  composite  laminate 
specimens  under  uniform  extension.  In  Section  II,  a  review  of  analytical  and 
numerical  methods  related  to  the  free-edge  stress  problem  is  presented.  Section  III 
contains  the  theoretical  foundation  of  the  finite  element  formulation  including  basic 
variational  principles.  Section  IV  describes  a  continuous  strain  finite  element  mcxlel 
based  on  a  compatible  cubic  interpolation  function.  A  continuous  traction  finite 
element  pnxedure  for  analysis  of  free-edge  delamination  specimens  is  developed  in 
Section  V.  In  Section  VI,  analysis  of  free-edge  effect  as  well  as  onset  of  delaminatinn 
in  various  types  of  laminated  specimens  are  presented.  Section  Vll  contains  discussion 
of  the  proposed  finite  element  models  for  analysis  of  free-edge  delamination  specimens. 
Derivation  of  Felippa's  compatible  cubic  interpolation  function  is  summarized  in  the 
Appendix. 
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SECTION  n 

REVIEW  OF  EARLIER  WORK 


2.1  Introduction 

'I'he  problem  ol'  caku bring  interlaminar  stresses  near  the  I'ree-edges  of  a  layered 
composite  under  uniform  inplane  extension  has  Iteen  investigated  by  many  researchers. 
Most  approximate  solutions  [2-7,10-19]  are  based  ujxtn  elasticity  theory  and  treat  the 
problem  as  a  generalized  plane  strain  case.  This  is  because  first  of  all,  the  classical 
and  even  many  of  the  refined  laminate  theories,  are  single-layer  theories  which  do  not 
account  for  local  effects  such  as  geometric  and  material  discontinuities,  and  the  presence 
of  a  free-edge;  secondly,  use  of  discrete  layered  theory  is  very  uneconomical  and 
impractical  from  the  computational  stand  point.  An  effective  modulus  formulation  [27] 
in  which  each  layer  is  characterized  as  a  homogeneous,  anisotropic  material  has  been 
widely  used  [1-19].  A  complex  state  of  stress  with  high  gradients  has  been  noticed 
[l9]  in  the  neighborhood  of  the  free-edge  due  to  the  presence  of  interlaminar  stresses  to 
keep  the  laminae  in  a  state  of  equilibrium.  In  order  to  have  a  precise  prediction  of 
deiamination  Itehavior,  an  accurate  estimate  for  the  near-field  stress  distribution  is 
essential.  However,  due  to  the  singular  nature  of  the  boundary-layer  stress  field  [19], 
an  exact  solution  is  currently  unavailable,  and  discrepancies  exist  in  the  magnitude  and 
even  the  sign  of  the  computed  interlaminar  stresses  near  the  free-edge  (Figure  1)  based 
on  various  approximate  theories. 
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2.2  Analytical  Approach 

Except  for  Pagano's  [8]  approximate  theory  based  on  Reissner's  variational  principle 
and  Pagano  and  Soni's  [lO]  Global-local  model,  most  analytical  solutions  discussed  in 
this  section  are  obtained  by  using  various  engineering  methods  to  solve  the 
displacement-equilibrium  equations  under  certain  assumptions.  Thus,  these  can  be 
regarded  as  approximate  solutions  based  upon  elasticity  theory. 

2.2.1  Approximate  Elasticity  Solution 

Investigations  ol'  the  I'ree-edge  problem  was  carried  out  by  Puppo  and  Evensen  [2] 
using  a  composite  model  es.sentially  consisting  of  a  set  of  anisotropic  layers  separated 
by  isotropic  adhesive  layers.  It  was  assumed  that  the  isotropic  layers,  developed  only 
interlaminar  shear  stresses,  acting  as  an  adhesive  between  the  anisotropic  layers.  It  w’as 
reported  that  a  sharp  rise  of  the  interlaminar  shear  stress  could  be  observed  in  finite 
width  laminates.  However,  the  simplicity  of  these  elastic  formulations  precluded 
calculation  of  the  transverse  normal  stress,  and  the  problem  became  more  complicated 
when  more  layers  were  involved. 

In  an  attempt  to  approximate  the  interlaminar  normal  stress,  a  simplified  formula 
was  developed  by  Pagano  and  Pipes  [l].  The  strategy  w'as  to  use  solutions  along  the 
longitudinal  mid -plane  of  the  laminate  based  upon  classical  laminated  plate  theory,  one 
could  then  compute  the  force  and  moment  re.sultants  caused  by  the  interlaminar  stresses 
on  any  plane  z=constant  through  consideration  of  static  equilibrium.  The  maximum 
interlaminar  normal  stress  at  the  free-edge  could  then  be  expressed  in  terms  of  the 
transverse  stress  in  the  y-direction  calculated  from  the  laminated  plate  theory. 

■Another  approximate  elasticity  solution  proposed  by  Pipes  and  Pagano  [3]  was  based 
upon  displacement-equilibrium  equations  for  an  anisotropic  elastic  medium.  Assuming 
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the  transverse  stresses  in  tlie  y-,  /-  directions  to  \anish,  the  equations  were  written  in 
terms  of  the  single  variable  U  (axial  displacement  function).  This  yielded  components 
of  displacement,  strain  as  well  as  remaining  stress  fields  in  the  form  of 
sinusoidal-hyperbolic  series.  However,  violation  of  stress  equilibrium  in  the  transverse 
directions  as  well  as  neglect  of  the  interlaminar  normal  stress  constituted  major 

drawbacks  of  this  scheme. 

2.2.2  Modified  Higher  Order  Theory 

Pagano  [4]  derived  another  approximate  metluxl  i'or  determination  ol'  distribution  of 
the  interlaminar  normal  stress  along  the  mid-plane  of  a  symmetric,  finite  width 
laminate.  The  approach  was  based  upon  a  modified  version  of  a  higher  order  theory 
proposed  by  Whitney  and  Sun  [28],  w'hich  recognized  the  effect  of  shear  deformation 
through  the  inplane  rotations  as  well  as  the  thickness  strain  implemented  in  the 
assumed  displacement  field.  However,  like  the  approximate  theories  discussed 

previously,  none  of  them  were  able  to  determine  the  complete  stress  field  near  the 
f  ree-edge. 

2.2.3  Galerkin's  Method 

Due  to  the  fact  that  high  stress  gradients  occurring  near  the  free-edge  are  difficult 
to  estimate  by  numerical  approaches,  Wang  and  Dickson  [5]  applied  the  extended 
Galerkin's  approach,  in  which  interlaminar  stre.sses  and  displacements  of  each  layer 
satisfying  geometrical  boundary  conditions  were  represented  as  series  of  Legendre 
polynomials.  The  final  solution  was  reached  by  requiring  the  satisfaction  of  continuity 

conditions  at  each  interface  as  well  as  stress  boundary  conditions  at  exterior  planes. 

Due  to  tl  '  completeness  of  I,egendre  polynomials,  convergence  of  solutions  could  be 
expected. 
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2.2.4  Perturbation  Technique 

In  an  effort  to  obtain  more  accurate  free-edge  stress  intensities,  a  perturbation 

technique  was  applied  by  Hsu  and  Herakovich  [6]  to  solve  the  three  coupled 

dimensionless  partial  differential  equations  based  upon  a  displacement  formulation  of 
the  elastic  problem.  It  showed  that  the  perturbation  solution  provided  a  smooth 
continuous  stress  distribution  in  the  vicinity  of  the  free-edge.  However,  this  solution 
had  the  limitation  that  the  shear  stress  distribution  was  a  function  of  both  laminate 

thickness-to-width  ratio  and  a  problem-de)>?ndent  parameter.  Although  the  latter  could 

be  chosen  sucli  that  the  maximum  values  of  shear  stress  field  did  not  exceed  elastic 
limits,  the  accuracy  of  the  calculated  stresses  was  suspect. 

2.2.5  Boundary  Layer  Theory 

A  boundary  layer  theory  for  laminated  composites  in  plane  stress  was  developed 
by  Tang  and  Levy  [7]  from  the  three-dimensional  theory  of  anisotropic  elasticity.  By 
expanding  the  stresses,  displacements,  body  forces  and  surface  tractions  in  power  series 
of  the  half-thickness  of  a  lamina  in  the  equations  of  equilibrium,  compatibility  and 
boundary  conditions,  a  sequence  of  systems  of  equations  was  obtained.  The  complete 
solution  was  obtained  by  combining  solutions  of  the  interior  domain  based  on  the 
classical  lamination  theory  and  those  from  Ixiundary  layer  and  matching  a  set  of 
appropriate  boundary  conditions.  This  formulation  indeed  provided  a  way  to  obtain 
analytical  solution  for  e.stimating  interlaminar  normal  as  well  as  shear  stress 
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2^.6  Reissner's  Variational  Principle 

In  order  to  have  displacement  as  well  as  stress  continuity,  a  mixed  formulation  is 
sometimes  used.  Unlike  the  elastic  approximations  discussed  previously,  Pagano  [8] 
developed  an  approximate  theory  for  a  general  composite  laminate  btised  upon  an 
application  of  Reissner's  variational  principle.  In  this  theory,  the  inplane  stresses  are 
considered  linear  in  the  thickness  coordinate  while  the  transverse  stresses  derived  from 
equilibrium  consideration  are  cubic.  If  a  laminate  or  a  single  lamina  is  viewed  as  an 
assembly  of  N  sheets,  each  having  a  finite  thickness  and  required  to  satisfy  force  and 

moment  equilibrium,  the  analysis  led  to  a  set  of  23N  algebraic  and  ordinary 

differential  equations  which  had  to  be  solved  simultaneously.  Based  upon  the 

assumption  that  the  stress  field  is  independent  of  the  longitudinal  axis,  Pagano  [9] 
further  specialized  the  theory  to  the  free-edge  problem  by  reducing  the  stress  field 
determination  to  the  solution  of  a  one-dimensional  problem.  Despite  the  relative 
accuracy  of  this  theory  resulting  from  the  improvement  of  smoothness  for  both 
displacement  and  traction  fields  at  interfaces  between  adjacent  layers,  a  major  drawback 
was  that  its  application  was  limited  at  most  to  six  sublayers. 

22..1  Global-local  Model 

Pagano  [lO]  introduced  a  global-local  model,  which  was  able  to  define  detailed 
response  functions  in  a  particular,  predetermined  region  of  interest  while  representing 
the  remainder  of  the  domain  by  effective  properties,  that  reduced  the  number  of 
variables  in  a  given  problem.  In  this  model,  for  the  global  region  of  the  laminate, 

potential  energy  has  been  utilized,  and  the  displacement  components  were  based  upon 

the  assumption  given  by  Whitney  and  Sun  [28].  The  Reissner  variational  principle 

described  in  [8l  however,  was  applied  for  the  local  region  in  which  a  thickness 
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distribution  0/  stress  I’ield  satisfying  equilibrium  equation  within  each  layer  was 
assumed.  A  variational  principle  was  then  used  to  derive  the  governing  equations  of 
equilibrium  for  the  whole  system.  It  was  reported  that  the  global-local  model  could 
effectively  solve  the  same  class  of  free-edge  stress  problem  as  described  in  [8]  and  had 
wider  range  of  applicability. 

2.3  Finite  Difference  Method 

2.3.1  Pesudo  Two-dimensional  Analysis 

Ihjres  and  I'agano  [ll]  used  the  classical  ilieort  ol  linear  elasticity  to  formulate  the 
problem  of  free-edge  delamination  ol  a  strip  under  unilorm  axial  strain.  .Allowing  for 
material  symmetries  and  uniform  extension,  the  transverse  components  of  displacement 
were  assumed  to  be  independent  of  the  longitudinal  coordinate.  The  three  coupled 
elliptic  equations  for  the  displacement  functions  were  solved  using  a  finite  difference 
solution  technique  to  approximate  the  interlaminar  stresses.  Delamination  was  assumed 
to  be  primarily  due  to  the  high  shear  stress  near  the  free-edge  and  the  interlaminar 
stress  field  was  found  to  be  an  edge  effect  which  was  restricted  to  a  boundary  region 
approximately  equal  to  the  laminate  thickness. 

2.3.2  Three-dimensional  Analysis 

A  three-dimensional  finite  difference  analysis  was  carried  out  by  Altu.s,  Rotem  and 
Shmueli  [12]  to  examine  the  free-edge  stress  field.  The  displacement  equilibrium 
equation  was  solved  by  using  central  difference  method  while  for  displacement  or 
traction-free  boundary  conditions  as  well  as  interfacial  continuity  conditions,  either 
forward  or  backward  difference  scheme  w'as  applied.  Convergence  of  the  solution  was 


n 


expecled  prDviding  a  reast)nable  displacement  lield  was  assumed  initially.  Although  a 
complete  stress  field  was  available  due  to  three-dimensional  characteristics,  an  iteration 
scheme  could  be  a  serious  inconvenience. 

2.4  Finite  Element  Method 

In  order  to  more  effectively  evaluate  the  high  gradient  stress  field  at  the  free-edge 
of  laminated  ci'inposites,  tlie  piipular  finite  element  method  has  been  applied  bv 
numerous  in\ estigaiors. 

2.4.1  Displacement  Method 

Wang  and  ('rossman  [13]  used  a  very  fine,  constant  strain  triangular  element  grid 
to  model  the  laminate  boundary  region  through  a  cross-section.  The  functional 
dependence  of  the  assumed  displacement  field  was  of  the  same  type  as  in  Pipes  and 
Pagano's  analysis  (ill.  To  overcome  the  difficulty  of  computational  storage  and  time 
limitation,  the  solution  process  adopted  the  so-called  ’’sky-line”  matrix  storage  scheme. 
The  results  indicated  that  the  interlaminar  as  well  as  inplane  stress  singular  behavior 
was  highly  localized  in  angle-ply  laminated  composite.  A  simplified  method  for 
calculating  interlaminar  stress  was  proposed  |14]  wherein  the  stresses  at  the  desired 
layer-interface  were  evaluated  by  substructuring  the  laminate  with  fewer  numliei  of 
effective  layers.  This  reduced  the  number  of  laminar  interfaces  and  facilitated  linite 
element  calculation  within  fewer  elements. 

A  qua.si-three-dimensional  finite  element  analysis  was  carried  out  by  Raju  and 
Crew  (15]  using  eight-noded  isoparametric  elements.  In  r>nler  to  approximate  the  stress 
singularities,  polar  mesh  was  introduced  near  the  intersection  of  interface  and  free-edge, 
associated  with  a  so-called  log-linear  procedure  to  relate  the  steep  gradient  stress  with 


12 


the  radial  distance  from  the  singular  point  in  the  logarithmic  coordinate.  One  major 
drawback  of  this  scheme  is  that  the  power  of  singularity  has  to  be  determined  by 
solutions  calculated  from  finer  polar  mesh  near  the  interface  of  the  free-edge. 

Whitcomb,  Raju  and  Goree  [16]  further  pointed  out  that  the  disagreement  for  both 
magnitude  and  sign  of  the  interlaminar  normal  stress  distribution  among  various 
numerical  methods  could  be  attributed  to  the  unsymmetric  stress  tensor  at  the 
singularity.  In  their  approach  too,  the  problem  was  modeled  by  eight-noded 
isoparametric  elements.  It  was  concluded  that  finite  element  displacement  models  were 
capable  of  giving  accurate  stress  distributions  everywhere  except  in  the  region  within 
two  elements  of  a  stress  singularity. 

In  summary,  we  observe  that  in  the  conventional  displacement-based  finite  element 
formulation,  evaluation  of  shear  as  well  as  normal  stresses  required  expensive  mesh 
refinement  near  the  boundary  region  to  approximate  the  singular  stress  field.  Even 
then  the  actual  stress  distribution  along  the  free-edge  was  generally  not  sufficiently 
accurate. 

2.4^  Stress  Method 

Rybicki  [l7]  used  a  three-dimensional  equilibrium  finite  element  analysis  procedure, 
based  upon  minimization  of  complemenury  energy,  to  solve  the  free-edge  stress 
problem.  Due  to  the  fact  that  the  assumed  stress  state  in  the  analysis  did  not  contain 
singular  term,  a  finite  rise  in  interlaminar  normal  and  shear  stresses  near  the  interface 
corner  was  observed  for  angle-ply  layup.  However,  this  method  involved  very  iarge 
matrices  and  was  computationally  expensive,  and  even  at  that  did  not  yield  a 
continuous  stress  field. 
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2.4.3  Hybrid  Assumed  Stress  Model 

In  J'ian's  hybrid  model  [29],  stress  equilibrium  in  the  interior  of  the  elements  as 
well  as  displacement  continuity  along  interelement  Ixtundaries  are  ensured,  but  the 
interelement  stress  continuity  is  satisfied  only  in  a  weighted  integral  sense.  Following 
Plan's  formulation,  Spilker  [IS]  developed  a  special  hybrid  element  for  the  edge-stress 
problem.  In  his  w'ork,  the  assumed  stress  field  was  made  to  satisfy  exactly  the 
continuity  ol  tr;u.tion  across  interlayer  Ixnindaries  as  well  as  traction-free  ronditit'O'^ 
along  exterior  planes  of  the  laminate.  This  was  I'ound  to  be  effective  for  study  of 
cross-ply  laminates  having  a  relatively  simple  stress  field.  It  is  difficult  to  extend  this 
procedure  to  angle-plv  laminates  because  in  these  the  complete  stress  field  has  to  be 
considered. 

A  special  formulation  of  a  .singular  composite-edge  element  was  developed  by 
Wang  and  Yuan  [l9]  based  on  the  Boundary-layer  theory  [30]  and  the  variational 
principle  of  a  modified  hybrid  functional.  In  the  analysis,  the  singular  hybrid  element 
was  used  in  conjunction  with  displacement-based  eight-noded  isoparametric  elements,  and 
it  w'as  reported  to  give  satisfactory  stress  distribution  near  the  free-edge.  This  method 
is  excellent  for  determining  possible  grow'th  of  delamination  but  would  be  awkward  to 
use  to  predict  CKCurrence  ol'  delaminalion  in  an  intact  specimen.  This  is  because 
sometimes,  it  is  hard  to  find  the  place  in  which  stress  singularity  may  (x;cur. 

2.5  Summary  and  Research  Motivation 

The  analytical  and  numerical  solutions  discussed  above  for  the  free-edge  stress 
problem  are  summarized  in  Table  (1).  Some  conclusions  can  be  made  at  this  point. 
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1.  Generally  speaking,  an  analytical  solution  lor  the  complete  stress  field  is 
extremely  difficult.  The  solution  pnxedure  for  the  case  of  a  multi-layer 
system  is  not  currently  available. 

2.  Use  of  a  finite  difference  technique  suffers  from  geometric  limitations. 

Calculation  of  stresses  at  the  interfaces  or  laminate  boundaries  needs  to  apply 
additional  techniques,  such  as  iteration  scheme.  Even  then  the  solution 

generally  lacks  edibility. 

,1.  Conventional  displacement-based  finite  element  methods  are  incapable  of 

predicting  accurate  stress  I  ields  particularly  along  element  boundaries.  Stress 
equilibrium  approach  is  apparently  impractical.  Use  of  hybrid  element  does 

improve  stress  calculation  but  is  applicable  only  to  some  s]tecial  cases. 

Application  of  singular  element  near  the  free-edge  boundary  apparently  makes 
the  analysis  too  subjective.  In  order  to  have  reliable  predictions  of 

displacement  and  stress  fields,  it  is  necessary  that  the  free-edge  stress  model  be 
able  to  approximate  the  real  life  situation  as  closely  as  possible.  In  other 
words,  the  displacement  and  stress  continuity  conditions  along  with 

traction-free  boundary  condition  have  to  be  exactly  satisfied.  Considering  also 
the  generality  and  effectiveness  of  the  analysis,  the  displacement-based  finite 

element  approach  with  higher  order  interpolation  function  could  conceivably  be 
superior  to  the  other  approximate  thcorie.s. 
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1'able  1:  Comparison  of  Various  Methods  for  Solving  the  I'ED  Problems 


Ref 

Author 

method  of  analysis 

calculated  stresses 

2 

Puppo  &  Evensen 

Elastic  approximation 

O’  a  T  T 

N  y  X2  xy 

3 

Pipes  &  Pagano 

Approximate  elastic  solution 

cr  7  T  r 

\  w  \7.  yv. 

4 

Pagano 

Modified  higher  order  theory 

a 

y 

5 

'A'ang  h  Dickson 

Extended  Galerkin's  approach 

a  T 

y  yy 

6 

Hsu  &  Herakovich 

Perturbation  technique 

(T  7  T 
z  \7.  yz 

7 

Tang  ti  Levy 

Boundary  Layer  theory 

O’  O’  a  T  7  7 

X  y  z  yz  X7  xy 

8 

Pagano 

Reissner's  variational 
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SECTION  III 

VARIATIONAL  FORMULATION  AND  FINITE 
ELEMENT  APPROXIMATION  IN  LINEAR  ELASTiaTY 

3.1  Introduction 

In  this  section,  a  variational  formulation  ol'  tliree-climensiona!  elasticity  is  tlescrilted 
and  its  use  as  the  basis  of  a  finite  element  approximation  is  discussed.  The  treatment 
es.sentially  follows  that  in  reference  [31].  Variational  formulation  has  been  used  as  the 
basis  lor  direct  methods  of  obtaining  approximate  solutions  to  boundary  value  and 
initial  boundary  value  problems.  Traditionally,  the  approximation  space  is  generated  by 
complete  orthonormal  sets  consisting  of  eigenfunctions  of  self-adjoint  operators.  The 
functions  which  are  used  to  approximate  the  field  variables  are  required  to  satisfy 
certain  continuity  requirements  over  the  whole  domain.  The  finite  element  method, 
however,  offers  an  alternative  route  for  generating  the  sequence  of  finite  dimensional 
approximation  spaces.  The  region  under  consideration  is  subdivided  into  a  finite 
number  of  discrete  elements,  and  the  field  variables  are  represented  by  functions  which 
follow'  the  stime  continuity  condition  only  piecewise  w'ithin  each  element.  Some 
significant  differences  between  the  finite  element  method  and  the  traditional  direct 
methods  include  [31] 

1.  The  base  functions  have  local  support  and  are  nonorthogonal. 

2.  The  sequence  of  approximation  spaces  is  ordered  by  refinement 

3.  The  l(x:al  support  functions  may  have  only  limited  smoothness 
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■J'he  support  of  each  base  function  is  confined  to  the  neighborhotxl  of  a  nixla) 
point  and  extends  over  the  elements  of  the  finite  element  approximation  sharing  that 
point.  Across  interelement  boundaries  within  the  support  and  at  the  boundary  of  the 
support,  the  function  may  have  only  limited  smoothness.  In  a  sequence  of  refinements, 
additional  nodal  points,  elements  and  base  functions  are  introduced.  The  base  functions 
ass(x:iated  with  each  ntxlal  point  change  with  refinement,  including  a  monotonic 
decrease  in  support.  Additional  discontinuities  might  Ire  introduced  at  each  refinement. 
Variational  formulations  and  solution  procedure  for  direct  methods  based  on  ilie  finite 
element  approach  must  allow  for  these  j^jculiarities  of  finite  element  approximation 
spaces. 


32  Boundary  Value  Problem 


Consider  an  open  connected  region  R  in  an  euclidean  space.  QR  is  the  boundary 
of  R  and  R  its  closure.  A  typical  boundary  value  problem  on  R  is  defined  by  the  set 
of  equations 


Au  =  f  on  R 

Cu  =  g  on  0R 

where  A  is  the  field  operator  and  C  is  the  boundary  operator  such  that 
A:D,-^V, 


(1) 

(2) 

'3) 

(4) 


Vr,  VgR  are  linear  vector  spaces  whose  elements  are  defined  on  the  regions  indicated 
by  the  subscripts.  D^,  Dg„  are  dense  subsets  in  V^,  and  denote  the  domains  of  A, 
C  respectively.  Dj„  is  the  extension  of  i.e.  any  element  u€Dr  has  a  unique 

extension  in  Dyp  and  every  element  in  is  the  extension  of  an  element  (not 
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necessarily  unique)  in  Dr.  For  given  fCVR,  g€VjR,  the  boundary  value  problem 
consists  of  determining  u€Dr  along  with  its  extension  in  Djr  such  that  (l)  and  (2)  are 
satisfied. 

33  A  Variational  Principle 

Let  the  linear  operator  A  be  self-adjoint,  i.e.  there  exists  a  nondegenerate,  linear 
Gateaux  differentiable,  bilinear  mapping  Br:DrXVr-*S,  where  S  is  a  linear  vector  space, 
such  that 

Br(u,Av)  =  Br(vAu)  +  C^r(v,u)  u,v€DgnVjj  (5) 

Here,  Cjr(v,u)  are  quantities  associated  with  the  boundary  QR.  Magri  [32]  has  shown 
that  such  a  bilinear  mapping  can  be  constructed  for  every  linear  operator  A.  If  the 
boundary  operator  C  is  consistent  [33]  with  the  field  operator  A,  i.e.,  there  exists  a 
nondegenerate,  linear  Gateaux  differentiable,  bilinear  mapping  B^:DjrXVjr-*S,  such  that 

C^(v,u)  =  B^(v,Cu)  -  B^(u,Cv)  (6) 

then,  the  linear  Gateaux  differential  of 

ft(u)  =  BR(uAu-2f)  +  B^(u,Cu-2g)  (7) 

vanishes  if  and  only  if  (l),  (2)  are  satisfied.  Sandhu  and  Salaam  [33]  further  pointed 
out  that  even  if  the  boundary  condition  is  homogeneous,  i.e.  g  =  0,  the  quantity 
Bjr(u,Cu)  in  (7)  must  be  included  if  the  variational  principle  is  to  hold  for  the  path 
of  Gateaux  differentiation  not  satisfying  homogeneous  boundary  conditions.  This  is 
important  for  approximation  in  finite  element  spaces  where  the  variation  is  introduced 
as  change  in  the  nodal  point  value  of  the  field  variable  and,  consequently,  the  path  of 
variation  may  not  satisfy  the  boundary  condition  or  internal  smoothness. 


3.4  Variational  Principle  for  Finite  Element  Approximation 


In  the  finite  element  method,  the  region  R  is  approximated  by  a  set  of  elements 
{R^;  e=l,2 . m)  such  that 

R^nRf  =  0  if  e;^f  (8) 

m 

lim  1)  R,  =  R  (9) 

III -*00 

t=l 

The  field  \ariables  are  appri'ximated  by  functions  w  hicli  may  not  be  sufficiently 
smooth.  However,  oxer  each  element,  ade(]uate  smcxithness  is  assured.  If  R^,  represents 
the  interior  of  the  e-th  element  and  QR,  its  Ixrundary,  we  have  [34] 

11^  (u,.Av)  =  (v,Au)  +  C  (v,u)  ( 10) 

and 

( v,u)  =  (v,Cu)  -  (u,Cv)  (11) 

Further  define 

m 

fi(u)  =  2^[Bj,  (u,Au— 2f)  +  (u,Cu— 2g)]  +  B  j(u,(Cu)')  (12) 

.=1  '  ' 

where  9R[  represents  interior  boundaries  of  elements  and  a  prime  denotes  a  jump.  The 
Gateaux  differential 

M) 

A^fi(u)  =  2  22[B^(v,Au-f)+B^^^^^(v,Cu-g)]  +  2B  ,(v,(Cu)')  (13) 

e=l 

vanishes  if  and  only  if  (l),  (2)  are  satisfied  over  each  element  and  (Cu)  vanishes,  i.e. 
Cu  is  continuous  across  interelement  boundaries.  If  there  are  actual  discontinuities  in 


the  interior  of  R,  let 

(Cu)'  =  g'  over  9^' 


(14) 


where  g'  is  specified  over  QR'.  Then,  if  the  union  of  intersection  of  element 
boundaries  covers  R',  the  functional  in  (12)  may  be  redefined  as 


n(u)  =  £lBg^(uAu-2f)  +  Bgg^^^(u.Cu-2g)]  +  B^,(u.(Cu)'-2g')  (15) 


3.5  Linear  Operator  with  Adjoint  Splitting 

Many  physical  problems  can  be  written  in  tlie  form  of 


Au  =  I'u  +  T  iri'u  =  f  on  R  (16) 

where 

F:W,-V,  (17) 

T:W^-X„  (18) 

E:X^-»Yj  (19) 

T*:Yj-Vg  (20) 

T*  is  the  adjoint  of  T,  i.e.  such  that 

Bj(u,Tv)  =  B^(v,T  u)  +  C^(v,u)  (21) 


Here  B^;W„xX^-*S  and  B^ : Wj,xV„-»S.  S  is  a  linear  vector  space  and  Br,  B^  are 
continuous  non-degenerate  bilinear  mappings.  C,  T  are  symmetric,  i.e. 


B^(u£v)  =  Bg(vJ-u)  (22) 

Bj(uJv)  =  Bg(vJ^u)  (23) 

Introducing  €,  <r  through  the  equations 

Tu-€  =  0  on  R  (24) 

E€  —  or  =  0  on  R  (25) 


(16)  can  be  written  as 
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Fu  + 1’  O'  =  f  on  R 

Combining  (24)  through  (26),  these  constitute  the  coupled  system 

F  0  T*]  u  f 

0  E  -1  €  =  0  on  R 
T  -1  0  O'  0 

If  the  inverse  of  E  exists,  let  G=E'‘.  Then,  combining  (24),  (25) 

Tu  —  Go'  =  0  on  R 
(26)  and  (28)  are  the  coupled  system 

F  T  u  f 

=  on  R 

T  -G  O'  0 

(29)  is  referred  to  as  the  complementary  form. 

For  an  operator  with  adjoint  splitting,  let  the  boundary  conditions  on  u,  o'  be 


C,u  =  g,  on  Sj  3R 

(30) 

(31) 

The  discontinuity  conditions  are 

(C,u)'  =  g'j  on  s', 

(32) 

(C20')'=g'2  on  (33) 

where  S',  and  S',  are  interior  surfaces  imbedded  in  the  intersection  of  finite  element 
boundaries.  C,,  Cj  consistent  with  T,  T*  implies  the  existence  of  bilinear  mapping 
such  that 

Bg  (u,Tv)  =  (v.t’u)  +  B  (v,C2u)  —  B  (u,C , v)  (34) 

*  e  S| 

where  S*  are  complementary  subsets  of  boundary  0R,  of  element  e. 
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The  function  governing  variational  formulation  of  (27)  for  finite  element 
approximation  is 

m 

n(u,€,o')  =  (uJ^u+T cr— 2f)  +  (eJEc— tr)  +  (tr, Tu—e)  +  B^^  (cr.CjU— 2g,)] 

e=i  '  *  '  *  ' 

+  [Bg„^^5^(u.C20-2gp]  +  B‘ (cr.(C,u)'-2g',)  +  jr)' -2g\)  (35) 

It  is  important  to  note  that  even  if  there  are  no  interior  discontinuities  in  the 
physical  problem  or  the  specified  boundary  conditions  are  homogeneous,  the  boundary 
and  the  discontinuity  terms  must  be  included  to  accommodate  the  nature  of  the  finite 
element  approximation  space  [32]. 


3.6  Principle  of  Minimum  Potential  Energy 

The  field  equations  for  isothermal  quasi-static  deformation  of  anisotropic,  linear 
elastic  solids,  assuming  no  initial  stresses  and  strains,  are: 

a)  Equilibrium  of  stresses 

o".^.  +  f^  =  0  on  R  (36) 

b)  Kinematics 

For  small  deformation,  the  strain-displacement  relationship  is 


c)  Constitutive  relations 


tr  .  =  E 

ij  ijkl  kl 


(38) 


on  an  open  bounded  connected  set  R  contained  in  the  three-dimensional  Euclidean  space 
E.  Here  u,,  f,,  €^,,  or,j,  E,jn,  are,  respectively,  the  components  of  the  displacement  vector, 
the  body  force  vector,  the  infinitesimal  strain  tensor,  the  symmetric  Cauchy  stress 
tensor  and  the  isothermal  elasticity  tensor.  The  range  of  indices  is  1,  2,  3  and 
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summation  on  repeated  indices  is  implied.  A  subscript  following  a  comma  denotes 
partial  differentiation  with  respect  to  the  coordinate,  in  the  reference  frame,  defined  by 
the  subscript.  Parantheses  around  subscripts  denote  the  symmetric  part  of  the  quantity. 

Let  the  functions  u,,  €,j,  ar^.  satisfy  the  continuity  and  differentiability  properties 
required  in  the  equations  of  elasticity  over  every  subregion  Rj.  Then,  admitting 
(u,,€jj,cr,j)  as  the  15- tuple  of  dependent  variables,  components  of  vectors  and  tensors 
being  regarded  as  ordered  subsets  in  an  n-tuple,  (36)-(38)  can  be  written  as  [33] 


0  0  il(s  A+s  A) 

2  aj  ai 

u 

1 

fk 

0  L.„,  -1 

ijkl 

= 

0 

i(s^  A+s  A)  -1  0 

2  ■'•ai  “ak 

(T 

1] 

0 

(39) 


Consistent  boundary  conditions  for  the  problem  are 


-n.u.  =  -n.fl. 

on  S, 

(40) 

<7.n.  =  t.  on 

U  >  J 

S2 

(41) 

where  the  nj  are  components  of  a  unit  normal  to  the  boundary  S,  and  the  jump 
conditions  are 


KV  =  g'2j  0“  s; 

-(n.u7  =  -g',..  on  S; 

Setting  up  the  problem  in  inner  product  space,  i.e. 


(42) 

(43) 


uvdR,  and  defining 

R 


B„(u,v)  =  2^Vu,v)j^ 

e=l 

The  basic  functional  corresponding  to  (35),  allowing  for  relaxed  continuity,  is  [33] 
n  ,(u,€„,o-..)  =  J  (€.£.^,6, -U^o-., -2€^^a. -2UT.+U, ^^.)dR 


(44) 
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f  r  u/(<7,jnp'-2g'2j)dS 

•»  Sj  *»  S,  J  s'^ 

-  /  (45) 

s; 


In  using  (45)  as  the  basis  for  finite  element  approximations,  it  is  not  necessary  for  the 
interpolants  to  satisfy  any  boundary  conditions  or  interelement  continuity.  For  no 
jump  discontinuities,  g',jj  and  g'j^  vanish.  However,  it  is  important  to  retain  the  terms 
containing  (c.^n,)',  (njU,)'  in  the  formulation. 

Using  symmetry  property  of  the  operator  matrix,  i.e. 


J^o-.^u.dR  =  -/^«r,.u,  dR  +  Ja.nn.dS  +  J  a.fn.u  US 

+  J 

to  eliminate  the  term  containing  cr,jj  from  (45),  the  functional  can  be  written  as 
n2(u.,€„,o-.)  =  J^€..E.^,€„dR  +  2  J^cr./u. -e.)dR  -  2  J^uf.dR 

—  2  r  uf  dS  — 2  r  <r.  n  (u  — fi  )dS  —  2  f  tr  . (n  u  )'dS 
Js,  '  J  s.  ^  ■  ■  J  si  ^  ■ 


(46) 


(47) 


is  the  modified  variational  principle  with  three  independent  fields  proposed  by 

Prager  [35].  If  Uj,  e,j  are  restricted  to  satisfy  the  last  of  (39),  the  strain-displacement 
relations,  Prager's  modified  principle  of  total  energy  theory  is  obtained 

n  =  f  €  E..,€,,dR-2  r  ufdR-2r  u.tdS-2r  o-  n(u-Q)dS 
3  J  „  U  OKI  kl  J  J  J  c  -J  J  ■  . 


—  if  a,  (nu  )'dS 

J  ‘J  J  ‘ 


(48) 
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(48)  \^as  also  proposed  by  Plan  and  Tong  [29]  and  is  the  basis  of  their  hybrid  method 
with  assumed  displacement  field.  If  the  displacement  field,  u,,  is  further  restricted  to 
satisfy  the  displacement  boundary  condition  (40)  on  S,,  fij  reduces  to 


=  //.'.“S-S J  0',/n,u,)-dS  (49) 

2  S, 

If  the  finite  element  interpolations  are  chosen  to  identically  satisfy  displacement 
continuity  across  S',,  the  last  term  in  vanished  and  (49)  becomes 


The  vanishing  of  the  variation  of  flj  with  respect  to  the  displacement  components  u, 
implies  the  satisfaction  of  equilibrium  equations  (36).  This  functional  corresponds  to 
the  classical  principle  of  minimum  potential  energy  which  is  customarily  used  as  the 
basis  of  the  finite  element  displacement  formulation  of  the  elastostatics  problems. 


3.7  Assumed  Displacement  Finite  Element  Formulation 

For  the  boundary  value  problem  stated  in  (l)  and  (2),  the  solutions  u  to  the 
forcing  functions  f  in  general  belong  to  the  space  of  square  integrable  function.  Lj 
is  a  separable  Hilbert  space.  However,  u  may  be  contained  in  a  subset  D  of  Lj  such 
that  A,  the  linear  operator,  is  defined  on  D.  We  assume  that  D  is  dense  in  Lj.  If 
the  set  of  functions  {0,^,  k=l,2, — oo}  is  a  basis  in  D,  then  any  function  uCL^  can  be 
expressed  as  an  infinite  sum: 

OO 

u  =  (51) 

k  =  l 
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A  scheme  to  generate  approximate  solutions  is  to  use  a  finite  set  of  terms  in  the 
infinite  sum  above.  Thus,  as  an  approximation 


u  =  (52) 

k=l 

The  approximation  process  consists  of  an  appropriate  choice  of  n,  (5^  and  the  coefficient 
a^.  Several  alternative  procedures  are  available.  The  finite  element  method  is  a 
special  process  of  selection  of  finite  subset  of  the  basis  {<^^}.  The  coefficients  a^  are 
generally  evaluated  by  requiring  the  approximate  solution  to  satisfy  a  variational 
principle. 

The  finite  element  idealization  essentially  partitiotrs  the  spatial  region  R  into  a 
finite  number  of  nontrival  discrete  elements  or  subregions.  The  geometry  of  the 
elements  is  defined  by  a  set  of  points  in  space  called  the  nodal  points  of  the  system. 
Over  an  element  m.  let  an  approximation  to  u  be  u"  such  that 


m  _  m. 


k  =  l 


(53) 


or  in  matrix  form,  dropping  the  subscript  n, 

u™  =  (54) 

where  is  a  row  vector  consisting  of  0™  as  its  elements  and  {a”}  is  a  column 

vector  of  coefficients  a”  Evaluating  the  function,  and  its  derivatives  up  to  a  certain 
order  at  nodal  points,  yields 

{un  =  [<f{a"*}  (55) 

where  {u”}  is  the  vector  of  nodal  point  values  of  the  function  and  its  derivatives  up 
to  the  order  selected,  and  is  the  matrix  of  base  functions  evaluated  at  each  nodal 
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point.  The  rows  and  columns  of  [^1”^  are  linearly  independent.  If  square,  the  matrix 
is  invertible.  Hence,  we  can  write 

{a”}  =  {u™}  =  [AF ‘  {u™}  (56) 

where  A=[^|"F 

Substituting  (56)  into  (54) 

u"*  =  ir f  [AT  ‘  {u;"}  =  I0"'f  {u™}  (57) 

where  [0"']  can  now  be  regarded  as  a  set  of  interpolation  functions  relating  nodal  point 
values  of  a  function  and  its  derivatives  up  to  a  preselected  order,  to  the  value  at  an 
arbitrary  point  within  the  element  m  [36]. 

In  applying  the  potential  energy  functional  shown  in  (50),  6,^  is  assumed  to  satisfy 
the  strain-displacement  relationship,  and  the  displacement  field  should  satisfy  the 
prescribed  displacement  boundary  condition  on  S,.  Vanishing  of  the  variation  would 
imply  satisfaction  of  the  equilibrium  equations.  As  stated  previously,  in  the  finite 
element  method,  the  displacement  u  is  approximated  by  interpolation  functions  and 
generalized  displacements  at  a  finite  number  of  nodal  points  of  each  element.  The 
interpolation  function  must  be  chosen  in  such  a  way  that  when  the  nodal  point 
displacements  for  two  adjacent  elements  are  compatible,  the  displacements  along  the 
common  boundary  are  compatible.  Meanwhile,  the  interpolation  function  must  also 
satisfy  the  requirement  that  the  first  derivatives  of  the  displacement  field  exist. 

Based  on  (57),  the  assumed  displacement  over  an  element  can  be  rewritten  in 
matrix  form  as 


(58) 
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where  q”  is  the  column  matrix  of  generalized  displacements  at  boundary  nodes  of 
element  m  which  determines  the  inter-element  continuity,  and  r”  is  the  column  matrix 
of  generalized  displacements  at  nodal  points  either  at  the  boundary  or  at  the  interior 
of  element  m  but  which  do  not  affect  the  interelement  continuity  requirements.  The 
corresponding  strain  distribution  is 


(59) 


where  <f>^  and  0"  are  obtained  by  differentiating  and  with  respect  to  the 
spatial  coordinates.  Substituting  (59)  into  (50)  and  expressing  in  the  finite  element 
discretized  form,  we  have  [37] 


T 

n 


where 


[£'"}=  matrix  of  elastic  constants  for  element  m 
[<^™}=  matrix  of  interpolation  functions  for  the 
body  forces  in  element  m 
[0™}=  matrix  of  interpolation  functions  for  the 

prescribed  tractions  on  the  surface  of  element  m 

The  summation  sign  in  (60)  implies  the  direct  stiffness  assembly  procedure,  and  the 
q 

vector  is  the  vector  of  global  displacements,  ftj  can  also  be  written  as 

M 

^  2  {q”}  +  {r”} 

m=l 

-2{F^"’f{q"’}-2{F;"}V})  (61) 
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where 


K  =  r  [0'"][E'"]t<^"’fdR 

qq  I  ‘•^eq-*  •*  ‘■^eq-*  m 

(62) 

K  =  r  [0’"][E'"][<^'“fdR 

rq  1  '"■^eq**  m 

m 

(63) 

K„  =  / 

m 

(64) 

F  =  r  [<^‘"][0n"{f"’)dR„+  r 

R™  ** 

m  /  m 

(65) 

F^=  r  [0:‘][0n^{f"'>dR„.+  f  [<][<Arf{t"’}ds„ 

R„  *'  s,ns„ 

fn  ^  ID 

(66) 

The  displacements  {r”}  in  element  m  are  independent  of  displacements. 

{r‘},  for  i^m. 

The  stationarity  condition  with  respect  to  their  variations  yields 

[K”]{q”}  +  [K;p{r"}-{F“}  =  0 

(67) 

Solving  (67)  for  {r“}  yields 

{r"}  =  [K;"j‘({F;"}-[K”]{q"‘}) 

(68) 

Substituting  (68)  into  (61)  yields 

M 

Oj  =  £({q'"f  [K™]  {q"’}  -  {F™f  {q”}  +  C„) 

m=l 

(69) 

where  [K™]  and  {F™}  are,  respectively,  the  element  stiffness  matrix  and 

the  equivalent 

nodal  forces  defined  by 

[K*"] = [K”  ]  -  [K;;f [k;"j  ‘[k^i 

(70) 

{f'"}={f”}-ik;"/[k;"j'{f;"} 

(71) 

=  -{F™}  [K”r*{F”}  =  constont 

(72) 
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(Ij  in  (69)  is  given  in  terms  of  the  generalized  displacements  {q}  which  are  not 
independent  for  different  elements.  Using  global  coordinates,  (69)  can  be  written  as 

n5  =  {qrtK]{q}-2{qf{F}  +  C'^  (73) 

Taking  the  variation  of  this  discretized  form  of  the  functional  yields  the  system  of 
algebraic  equations 

[K]{q}  =  {F}  (74) 

which  can  be  solved  for  the  unknown  nodal  displacement  {q}.  The  matrix  [K]  is 
positive  definite,  symmetric  and  banded.  The  process  of  eliminating  the  generalized 
coordinates,  {r},  from  each  element  is  called  the  static-condensation  process  [38].  The 
introduction  of  these  terms,  which  do  not  interfere  with  interelement  compatibility, 
results  in  an  improvement  in  the  satisfaction  of  the  equilibrium  equations  within  each 
element.  However,  the  satisfaction  of  the  equilibrium  equations  along  the  interelement 
boundary  is  still  governed  by  the  degrees  of  compatibility  supplied  by  the  interpolation 
functions  for  the  generalized  displacements,  {q}.  The  solution  obtained  represents  an 
underestimate  of  the  true  solution  in  the  sense  of  energy  [39]. 
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SECTION  IV 


CONTINUOUS  STRAIN  FINITE  ELEMENT 
INTERPOLATION 


A.\  Introduction 

In  the  finite  element  metluKl,  tlie  displacement  lield  is  approximated  l)v 
interpolation  functions  and  generalized  displacements  at  a  finte  number  of  ntxlal  jxiints 
which  also  del’ine  the  geometry  of  the  elements,  'lb  ensure  continuous  strain  across 
interelement  boundaries,  it  is  sufficient  that  the  interpolation  functions  be  such  that 
the  displacement  components  as  well  as  their  first  derivatives  along  the  common 
boundary  are  continuous. 

Tocher  and  Hartz  [40]  pointed  out  that  for  plate  bending  analysis,  continuity  of 
slopes  of  the  plate  displacement  surface  is  necessary.  The  compatible  cubic  interpolation 
functions  developed  by  Tocher  [40]  and  by  Clough  and  Felippa  [4l],  among  others, 
satisfy  this  requirement.  For  plate  bending,  the  generalized  displacements  used  were  w, 
the  transverse  displacement  of  the  plate  and  its  derivatives  w,.,  w^.  Here,  the 
subscripts  x  or  y  denote  partial  differentiation  with  respect  to  the  independent 
variables  x,  y.  Applying  the  same  displacement  interpolation  scheme  to  the  plane 
elasticity  problems  [40],  the  corresponding  generalized  displacements  were 
u,  u„,  Uy,  V,  V,,  Vy,  the  in-plane  displacements  and  their  first  derivatives  at  each  node. 
Thus,  continuity  of  strain  between  adjacent  elements  was  ensured. 
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Instead  of  using  a  local  ('artesian  ciwrdinate  system  in  the  derivation  of  the  cubic- 
polynomial  with  nine  coefficients  for  each  displacement  component  in  Tocher  and 
Hartz's  -work,  triangular  coordinates  were  used  in  the  present  fully  compatible 
quadrilateral  element,  following  Felippa's  [42]  work  on  plate  bending  analysis.  This 
simplified  the  generation  of  various  matrix  relationships  for  the  constituent  triangular 
elements.  Tocher's  [40]  element  used  incomplete  cubic  polynomials.  Felippa's  elements 
-vvere  based  on  complete  cubic  poiviK'inials  .-incl  were  therefore  selected  for  application  to 
the  free-edge  pnrblems.  Tliese  elements  include  focher's  formulation  as  sj-tecialization. 
('ubic  e.xpansion  ol  the  9-degree-ol -freedvtm  conl'urming  triangular  element  (LCCr-9)  for 
lx)th  in-plane  displacement  components  as  used  by  Tocher  [40]  was  extended  to 
quadrilateral  element  designated  (}-l5.  (Quadrilateral  elements,  Q-19  and  (Q-23, 
assembled  from  LCCT-11  and  LC(Jr-12  triangular  elements  introduced  by  Felippa  [42], 
were  also  redeveloped  for  the  plane  elasticity  problems.  The  continuous  strain  elements 
were  used  to  analyze  a  pseudo  two-dimensional  free-edge  stress  problem  similar  to  that 
of  Pipes  and  Pagano  [II]  for  composite  laminate  coupons  under  uniform  extension. 

4,2  Interpolation  Functions  of  Continuous  Strain  Elements 

In  the  following.  Felippa's  [41,42]  approach  for  deriving  the  plate  bending 
interpolation  functions  is  summarized.  We  use  the  same  element  name  as  Felippa's  and 
start  with  u  in.stead  of  w  for  the  plane  elasticity  problems.  Similar  derivations 
applied  to  the  displacement  component  v. 
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4^.1  LCCT-12  Element 


A  complete  cubic  polynomial  m  two  variables  is  defined  by  ten  independent 
coefficients.  I'he  values  of  u,  the  x-direction  displacement  of  the  plane  stress  body  and 
its  derivatives  u^.,  u^,  at  the  three  vertices  of  a  triangle  yields  nine  independent 
quantities.  To  ensure  continuity  of  derivative  u„  across  element  boundaries,  it  is 
neces.sary  that  u„  be  known  at  some  points  other  than  the  vertices  along  each  of  the 
three  edges.  )i  is  consenient  to  introduce  mid  side  nodes  on  each  of  the  three  edges. 
This  element  then  has  twelve  indeisendent  quantities  against  the  minimum  of  ten 
needed  to  completely  define  a  cubic  polynomial. 

In  order  to  use  a  cubic  polynomial  with  continuous  first  derivatives  in  the  interior 
as  well  as  on  the  element  boundaries,  I'elippa  proposed  that  the  element  be  made  up  of 
three  subtriangles  as  illustrated  in  Figure  (2).  Fach  subtriangle  has  three  vertices  and 
one  mid-side  node  to  supply  the  ten  independent  quantities  for  defining  the  cubic 
polynomial  interpolation  in  its  interior.  The  point  O  could  be  any  interior  point. 
However,  for  simplicity  of  formulation,  the  centroid  is  generally  used  [41]. 

The  nodal  displacement  degrees  of  freedom  to  be  considered  in  the  stiffness  matrix 
of  the  complete  element  (Figure  2)  include  the  values  of  the  in-plane  displacement 
components,  u,,  v,  along  with  their  first  derivatives  u^^,  u^.,  v^^,  (i=l,2,.1)  about  the 

X  and  y  axes  at  each  corner  as  w^ell  as  the  normal  slopes  at  the  three  mid-side  ncxles 
about  axes  perpendicular  to  these  sides  respectively,  viz.  u,^,  u^^,  u^^,  and  v„^,  v^^,  v^^. 

After  forming  the  expression  of  the  cubic  displacement  patterns  in  the  three 
subelements,  because  of  the  common  displacements  imposed  at  the  nodes,  the  in-plane 
displacements  of  two  adjacent  subelements  are  identical  along  their  juncture  line.  To 
establish  continuity  of  u^  along  the  edges  of  the  subelements,  it  is  sufficient  that  u„ 
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evaluated  at  joints  7,  8,  9,  mid')xanls  ol  these  edj>es,  from  adjacent  subtriangles  be  the 

same.  These  three  conditions  were  used  to  evaluate  the  in-plane  displacement  u,  and 

its  derivatives  u^,  u^  at  the  interior  point  O.  With  the  interior  point  thus  condensed 
out,  Felippa  [41]  obtained  a  set  of  interpolation  functions  for  the  LCCT-12  element. 
These  define  a  piecewise  cubic  polynomial  interpolation  such  that  the  in-plane 
displacements  and  their  first  derivatives  are  continuous  both  in  the  interior  of  the 

elenieni  and  along  the  entire  boundarv  of  the  complete  triangular  element.  A  more 
detailed  deri\ation  procedure  and  the  complete  listing  of  the  cubic  interpolation 
functions  are  given  in  Appendix. 

4.2.2  LCCT-11  and  LCCT-9  Elements 

Assembly  of  three  subtriangles  results  in  the  LCCr-12  element  (Figure  3a). 
However,  the  mid-side  nodal  points  in  this  element  are  not  desirable  for  programming. 
They  complicate  the  mesh  generation  procedure,  increase  the  band-width  of  the 

assembled  equation  systems,  and  require  special  identification  in  calculation  of  the 
stiffness  matrix.  To  overcome  these  difficulties,  it  may  be  desirable  to  develop  a 
special  element  without  external  midpoints.  This  can  be  accomplished  by  assuming  the 
normal  slope  to  vary  linearly  along  one  or  more  sides  [42]. 

With  the  elimination  of  one  mid-side  node,  the  five-node  element  is  designated  as 
1,(XT-11  (Figure  3b).  Further  imposing  linear  slope  variation  constraints  on  three  sides 
gives  a  triangle  with  three  nodal  points  and  results  in  LCCT-9  element  as  illustrated 
in  Figure  3(c).  The  LCCT-9  element  is  identical  to  the  Tocher  and  Hartz  [40]  element. 
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4.2.3  Quadrilateral  Elements 


lilements  of  quadrilateral  shape  can  be  set  up  as  assemblages  of  triangular 
elements.  I'igure  (4)  shows  quadrilateral  elements  built  up  from  four  l.CCf-12, 
LCCT-ll  and  LCCT-9  elements.  The  quadrilateral  element  in  Figure  4(a)  has  a  total 
of  23  degrees  of  freedom  for  each  variable  and  was  designated  by  Felippa  as  Q-23. 
Using  four  l-Ud’-ll  or  l,CCl'-9  triangles,  the  (^-19  and  Q-15  elements  as  shown  in 
Figures  4(15).  4(t)  respetri\ elv  are  reali/.ed. 

The  (juadriluieral  element  lias  interior  nodal  points  not  connected  to  the  other 
quadrilateral  element  in  a  finite  element  mesh.  These  jxiints  can  be  eliminated  through 
a  local  condensation  prtKess.  Thus,  the  final  quadrilateral  element  has  24  degrees  ol 
freedom,  corresponding  to  the  two  in-plane  displacement  components  and  their  first 
derivatives  with  respect  to  spatial  coordinates  x  and  y  at  the  four  corners  of  the 
elements  and  an  additional  eight  degrees  of  freedom  corresponding  to  the  normal 
derivatives  of  each  of  the  displacement  components  at  the  mid-side  nodes.  A  uming 
that  the  normal  derivatives  vary  linearly  along  the  edges  of  the  quadrilateral,  the 
mid-side  nodes  can  be  dropped.  This  reduces  the  0-23  to  Felippa's  (J-19  element  with 
12  degrees  of  freedom  for  each  of  the  displacement  components.  It  is  a  fully 
compsilible  quadrilateral  element,  having  a  continuous  cubic  variation  of  displacement 
and  quadratic  variation  of  strain  both  in  the  interior  of  the  element  and  along  the 
entire  boundary  of  the  element,  as  well  as  a  linear  variation  of  normal  slope  along  all 
external  edges.  We  note  however  that  LCCT-9  and  LCCT-ll  do  not  use  a  complete 
cubic  polynomial.  For  this  reason,  Felippa’s  LCCT-12  element  based  on  complete  cubic 
interpolation  was  considered  an  improvement  upon  Tocher's  [40]  LCCT-9. 


38 


(a)  0-23  (b)  0-19 


(c)  0-15 


Figure  4: 


Ouadrilateral  Elements  formed  from  (a)  LCCT-12  (b)  LCCT-11  (c) 
LCCr-9 


4.3  Application  to  the  Free-Edge  Stress  Problem 

I'igure  (5)  shows  a  symmetric  laminated  com|x>site  coupon  under  a  state  of 
uniform  axial  strain.  In  this  case,  away  from  the  ends,  the  transverse  x=constant 
plane  displacement  fields  can  be  assumed  to  be  independent  of  x.  These  assumptions 
imply  the  following  form  for  the  three  components  of  displacement  [ll]. 
u(x,y,z)  =  ex  +  U(v,7.) 

v(x,y,z)  =  \’(y./)  (75) 

\\(.\,v,z)  =  \V(y,/) 

where  e„  is  the  unilorm  in  plane  strain  in  the  x-direciion  and  u,  v,  w  are  components 
of  displacement  along  x,  y  and  /  axes  res|x‘cti\ ely. 
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4.3.1  Finite  Element  Formulation 

The  constitutive  relationship  for  linear  elastic  anisotropic  material  obeys  the 
generalized  Hooke's  law 

o-=Ce.  i,j=l,2,....6  (76) 

where  e,  is  namely  the  uniform  extensional  strain  e„.  Based  upon  the  minimum 
potential  energy  principle  (50)  and  substituting  various  interpolation  functions  for 

displacement,  Ixidy  force  and  traction  fields  ap]'>earing  in  the  governing  functional  (60). 
a  nodal  rorce-displacemenl  relation  within  each  element  is  expressed  as 

K  u  =11  (77) 

l.l  I  I 

where  R,  represents  the  resultant  external  ntxial  force,  is  the  element  stiffness 
matrix  which  can  be  written  as 

K  =  r  B  C  B  dV  m,n=l,2 _ 6  (78) 

IJ  J  ^  iin  mn  nj  ’  ’  ’ 

The  range  of  i,  j  deytends  upon  the  'degree  of  freedom'  of  the  element,  B  is  the 
displacement  transformation  matrix  and  V  is  the  domain  of  the  element. 

Because  of  the  longitudinal  extensional  strain  is  specified  as  constant,  the 

corresponding  term  in  the  stiffness  matrix  can  be  separated  from  the  rest  and  (78) 
rewritten  as 

K  u  =  K  -  (79) 

M  J  >  • 

w’here  the  range  of  summation  on  m,  n  is  now  2,  3,  6  and  R°  is  the  element 

residual  force  due  to  uniform  in-plane  strain  e„,  i.e. 


R  =  /  B  C  ,edV 

i  I  im  ml  o 
V 
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After  forming  the  system  stiffness  matrix  and  nodal  force  vectors,  the  displacement 
components  can  Ire  obtained  by  solving  the  resulting  set  of  linear  equations  in  the 
standard  manner. 

4.32  Higher  Order  Elements 

For  the  free-edge  stress  problem,  due  to  the  fact  that  dependence  of  the 
longitudinal  displacement  on  the  longitudinal  ccxrrdinate  x  is  made  explicit,  the  three 
displacement  components  are  completely  delined  by  three  functions  of  two  independent 
trans\ersc  coordinates  y  and  z  as  shown  in  (75).  Thus,  the  compatible  cubic 
interpolation  functions  used  I'or  plane  elasticity  problems  can  Ire  extended  to  the  pseudo 
two-dimensional  model  of  a  laminate  couiwn,  and  a  continuous  strain  field  along  both 
in-plane  and  transverse  directions  ensured. 

Figure  (6)  shows  the  nodal  displacement  degrees  of  freedom  considered  in  the 
stiffness  matrix  for  the  complete  triangular  element.  These  included  the  values  of  the 
in-plane  displacement  components  u,  v,  the  transverse  displacement  w,  along  with  the 
first  deriviatives  u^,  u^,  Vy,  v^,  w^  about  the  y  and  z  axes  at  each  corners  i=l,2,3 
as  well  as  the  normal  slopes  at  the  three  mid-side  nodes,  viz. 

Unv  “n;,.  Un,-  v„^,  v^^,  w^^,  w„^.  Tliis  is  the  LCCT-12  element  but  with  total 

of  .16  degrees  of  freedom.  Further  assuming  the  normal  slope  to  vary  linearly  along 
one  or  all  three  sides,  the  LCXrr  il  and  LCXTl'-b  elements,  with  total  of  33  and  27 
degrees  of  freedom  respectively,  are  obtained  as  specializations. 

As  described  for  plane  elasticity,  quadrilateral  elements,  designated  as  Q-23,  Q-19 
and  Q-15,  were  set  up  as  assemblage  of  four  LCCT-12,  LCCT-ll  and  LC(7r-9 

respectively.  After  eliminating  the  interior  nodal  points  through  a  local  condensation 
process,  the  final  quadrilateral  element,  (^-23,  had  36  degrees  of  freedom,  corresponding 


to  the  three  displacement  components  and  their  first  derivatives  with  respect  to  the 
spatial  ccx)rdinaies  y  and  l  at  the  four  corners  of  the  element  and  an  additional  12 
degrees  of  freedom  corresponding  to  the  normal  derivative  of  each  of  the  displacement 
components  at  the  mid-side  nodes.  It  is  a  fully  compatible  quadrilateral  element, 
having  a  continuous  cubic  variation  of  displacement  and  quadratic  variation  of  strain 
not  only  within  the  elements  as  well  as  along  element  boundaries,  but  also  across 
laminate  interlaces. 
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SECTION  V 


CONTINUOUS  TRACTION  FINITE  ELEMENT 
PROCEDURE  FOR  COMPOSITE  LAMINATES 

5.1  Introduction 

In  the  continuous  stniin  Q-23  element  developed  for  the  free-edge  stress  problem, 
both  displacement  and  strain  are  continuous  along  interelement  as  well  as  interlaminar 
Ixiundaries.  However,  the  tractions  calculated  across  the  interfaces  between  differently 
oriented  layers  are  discontinuous  due  to  different  orientation  of  adjacent  plies.  Also, 
traction-free  boundary  conditions  associated  with  the  finite-width  laminate  coupon 
cannot  be  satisfied.  In  order  to  remedy  these  two  defects  and  make  the  numerical 
model  more  representative  of  the  real  situation,  it  was  necessary  to  ensure  interelement 
as  well  as  interlaminar  continuities  of  tractions  and  the  traction-free  boundary 
condition  along  the  free-edge.  To  accomplish  this,  nodal  point  degrees  of  freedom  must 
include  some  components  of  stress  and  exclude  normal  gradients  of  displacement  which 
will  Ire  different  across  interelement  Ixiundaries.  This  was  implemented  by 
transforming  the  displacements  and  their  normal  gradients  at  each  of  the  nodal  points 
of  the  Q-23  element  to  a  mixed  set  of  degrees  of  freedom  which  would  be  continuous 
across  interelement  boundaries.  These  included  both  displacement  and  interlaminar 
traction  components.  Appropriate  displacement-stress  relationships  derived  from  the 
constitutive  laws  were  used.  For  this  element,  traction-free  boundary  condition  could 
be  satisfied  in  a  point-wise  sense. 
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'I’he  continuous  traction  (J-23  element  still  lias  cubic  variation  of  displacement  over 
the  element  and  retains  continuity  of  displacements  across  inlerelement  as  well  as 
interlaminar  boundaries.  The  strains  as  well  as  stresses  vary  quadratically  within  each 
subtriangle  of  the  constituent  LCCT-12  elements  of  the  quadrilateral.  However,  certain 
components  of  strain  are  not  continuous  across  interelement  boundaries  but  interelement 
tractions  are.  This  correctly  allows  for  possible  differences  in  orientation  of  adjacent 
laminae.  II'  adjacent  layers  have  the  same  stress-strain  relationshi]-)s  due  to  identical 
orientation,  stress  continuity  will  imply  strain  continuity  as  well. 


5.2  Derivation  of  Displacement-Stress  Transformation 

In  the  Q-23  element  analysis,  the  number  of  nodal  degrees  of  freedom  is  different 
for  the  corner  nodal  points  and  the  midside  nodes.  For  this  reason,  derivation  of 
transformation  matrices  for  displacements  and  their  gradients  at  the  corners  of  the 
LCCT-12  element,  and  for  normal  ;radients  of  displacement  at  the  midside  nodes  on  the 
element  boundaries,  is  discussed  separately  in  the  following  sections. 

5.2.1  Corner  Nodes 

The  strain-stress  relationship  for  an  orthotropic  material  expressed  in  the  x-y-r. 
coordinate  system  is 

e  =  S  <T  i,j=l,2,....6  (81) 

or  in  matrix  form 
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where  S,j  are  components  of  the  compliance  matrix  for  monoclinic  materials  which 
have  symmeirv  with  respect  to  x-v  plane  (I'igure  2)  and  are  defined  as  [4't] 


S,  ,=S,  ,m''+(2S,,+S,_,  )m'n‘+S,2n'^ 

3  =  ^1. "■'■^23"^ 

Sj^=[2S^jm^-2S^2’’^‘'‘^2Sj2+S^^Xn^-m^)]mn 


S22=S,,n'*+(2Sj2+Sj^(_)m^n^+S22in'* 

523=S,3"'+S23m' 

§2^=[2S,,n^-2S22m4(2S,2+S(,6Xm^-n^)]mn 

533=^33 

S36=2(S,3-S23)mn 

S44='%4"^^+S53n" 

S.,=— S_mn+S,^mn 

45  44  55 

Ss5=S,,n'+S,^m' 

, ,  +S2  ,-2S ,  )m  "n 


(83) 


w'here  m=cos0',  n=sin0',  and 


^  V  V 

<  =_L.  s  =-- _1L  S  =-_iL 

^11  p  ’‘^12  p  ’^13  p 

*'11  *'22  *'33 


Q  ___  12  c  —  ^  c  — _  3  2 

*21~  p  ’  ^22~  p  *  ^23~  p 

'11  22  33 


(84) 


c  C  -  ^23  c  -  1 

■^31  p  *  ^32  p  *  ^33  p 

^11  ^22  ^33 


c  ^  c  —  ^  c  —  ^ 

^44  G  ’  ^55  G  ’  G 

''*23  ^13  12 

6'  is  angle  of  ply  from  global  axis  x  to  material  axis  1,  and  I-^^,  G,j,  (i,j,k=],2..T) 
are  moduli  of  elasticity,  shear  moduli,  l‘oisson's  ratios,  respectively,  in  material 
coordinates. 

Replacing  the  strain  components  by  their  corresponding  displacement  gradients,  for 
small  strain  theory,  and  rearranging  the  constitutive  relation,  (82)  becomes 
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or  symbolically 
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(85) 


(86) 


where 


k.}= 


{tr,}= 


u 
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II 

w  +v 
y  z 

u 

u 

y 

Z 

{(T,}  = 


.VZ 


(87) 


49 


and 


^1  1 
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0 

S4.S 
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(88) 


To  relate  the  interlaminar  strain  components  {e^}  to  tlieir  corresponding  interlaminar 
stress  components  {u,},  (a,}  was  eliminated  through  a  static  condensation  process.  'Ihis 
yields 

<el=fS]<a,}  (89) 

where 

{€}={€2HD2,]  [D,, ]  '(€,}  (90) 

[S'HD2,HD2,]  [D,,r‘[D,2]  (91) 

The  inverse  of  [D„],  namely,  the  compliance  matrix  in  plane  stress  case,  can  be  written 
explicitly  as 


Id,,) 


On  OnQ,j 

Q|2  ^22  ^26 

^26  ^h6 


(92) 


where 
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(93) 


^1 1=^1  ,m'‘+2(Q,2+2Qf,h)m‘n^+Q22n'‘ 

0,<,=-mn^Q22+m^nQ,j-mn(m^-n^XQ,2+2Q^^) 
^22=0.  ^nU2iQ^,+2Qjm^n^+Q^y 
^2^=— m^nQjj+iTin^Qj ,  +mn(m^— n^XQ,  2'*'2Q^^) 

066=(Qi.+Q22-2C»,2>"'"'+06> 

with 


\-v 


1 1 


21*^12 


-.  0„= 


'22 


1  ~V  V 
12  21 


Q,.=  .  o,,=(;.2 


*  2  ]  p 

12  21 


Substituting  (92)  into  (91)  and  (90),  (89)  could  be  expressed  as 


w  — B,u  — B,v  — B,u 

S,-X  0  0 

or 

z  lx  2  y  3  y 

33 

z 

W  +V 

2 

T 

y  z 

A4  45 

yz 

u 

0  5,^ 

T 

z 

45  55 

xz 

w’here 


16 


and 


B,=5„Q,.+5j,Q„+Sh<J 

®2“^1  3^  1 2'*'^23^22'*’^3<P26 
®3“^  1  3^  1 6  "^^2  3^2(.  '^^3(M,6 


X=BiS.3+B2S,3+B3S3, 


(94) 


(95) 


(96) 

(97) 

(98) 

(99) 


The  gradients  of  displacement  appearing  in  the  expression  for  interlaminar  strains  were 
then  written  in  terms  of  interlaminar  stresses  using  (95),  i.e. 

(100) 
(101) 


U  =S.,T  +S..T 
z  45  yz  55  xz 


V  =(w  +V  )— W  =S,,T  +S.,T  — W 

z  V  7  y  44  vz  45  xz  y 
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Combining  (1(X))-(102)  with  the  rest  of  the  displacement  nodal  degrees  of  freedom,  u, 
V,  w,  Uy,  Vy,  Wy,  and  noting  that  u^=e„,  the  applied  strain  loading,  the  generalized 
displacement  components  at  the  corner  nodes  of  the  LCC1’-12  element  were  related  to  a 
mixed  set  of  degrees  oi'  freedom  as  follows: 
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or  symbolically 

{r}=[G]{r'}+{R}  (104) 

B,,  Bj,  Bj,  and  X  occurring  in  (103)  have  been  defined  by  (96)  through  (99).  Thus,  at 
each  of  the  three  corners  in  the  I,C(7f  12  element,  we  have  three  displacement 
components  u,  v,  w  along  with  three  inplane  strain  components,  u^,  Vy,  Wy,  and  three 
interlaminar  stress  components  Ty^,  cr^. 


SJ22,  Mid-side  Nodes 
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three  traction  components  which  indeed  can  be  regarded  us  the  out-ol -plane  shearing 
stress  cr„^,  the  normal  stress  and  the  in-plane  shearing  stress  cr„..  'I'lie  relation 

between  Xj  and  x\  coordinates  as  shown  in  Figure  (7)  yields 


n|=0.  n,=m,  n^=n  (111) 

where  m=cos6,  n=sin0  and  9  is  the  angle  between  these  two  coordinates.  Substituting 

the  alx)ve  quantities  into  (109),  we  have 

^'2=°'nn='"\+2mnTy^+n"cr^  (113) 

t’3=(r^^=(m^-n^)T^^+mn(tr^— o-p  (114) 

Here,  for  interlaminar  stresses  acting  on  the  element  boundary  are  nothing  but  the 
traction  field,  <r„^,  tr„„,  cr„,  expressed  in  (112)-(114)  indeed  represent  the  three  traction 
components  as  illustrated  in  Figure  (7). 

The  stress-strain  relationship  for  the  orthotropic  lamina  expressed  in  the  x-y-z 
laminate  coordinate  system  is 
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(115) 
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where  is  tlie  stillness  matrix  for  monoclinic  symmetry  with  respect  to  x-y  plane 
and  is  defined  as  [43] 


Cj  |=C,  ,m"'+2(C,  ,+2C^^)m'n’+C,,n^ 


66 
2  "> 


Cj2=(C,,+C22~4C^^)m‘n‘'+C|2(w>  +n  ) 

Ci3=C,3m4C23n" 

C_,3=C,,n''+2(C,,+2q^)m^n-+C32m'' 


(',,11111  +('||ni  n— (('i  ,+  2C^  ^  'ninlm  -n  ) 

(:,,=-(:,,mSi+(; 


(116) 


(r,j  =—(),, m  •i+('|  ,11111  +((’,  ,+^(  )nin(m'— 11 


^46=(C,,s-<:44)'^'^ 
^55“^S5^  "^^44*^ 


C,  =(C , ,  +C , ,-2C ,  ,)m “"n^+C , (m ^-n 

00  11  c£  12  oO 

where,  with  m=cos0',  n=sin0',  again  0'  is  angle  between  the  global  axis  x  to  material 
axis  1,  and 

(l — p  )F'  (p  p  )F 

"11  A  ’  '12  ^ 


t^l,= 


with 
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.(’-^2^'2, 
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=G  ,, 

23’  5.S 

13  ’ 

— X/  V  — 1/ 

p  — ; 

12  21 

23  32 

'  1 1 


2.1 


(117) 


(118) 


(112)-(114)  along  with  (115),  gives 
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m^C,(_+n‘C3,^  m‘C2,+n’C33 

2mnC . . 

44 

(Cj^— C2g)nin  (C23— (m^— n'X! 


44 


2mnC^, 

4y 


2mnC'  , 

44 

(m^— n’)C 

44 


mC36  ' 

m^C,3+n-C33 

(C^^— r3,)mn 


+ 


mC,,e 

l(i  o 


j(niT  ,+n"C  )e 

I  1  -  I  ''I 

((', ('1  ,)mne^ 


(119) 


(119)  is  the  general  relationship  between  traction  components  at  the  element  niiti-sitie 
nodes  and  the  corres[X)nding  displacement  gradients.  It  is  difierent  for  each  element 
midside  node. 

In  order  to  relate  the  displacement  gradients  at  the  mid-point  to  the  mixed  nodal 
degrees  of  freedom  at  the  two  ends  defining  the  element  surface,  let  j  and  k  denote 
the  first  and  second  cyclic  permutations  of  i=l,2,3  (i.e.  j=2,3,l  and  k=3,l,2),  the 
projected  dimensions  and  the  corresponding  boundary  length  are  defined  as  (see 
Appendix  A) 


a  -y—y  ,  b  =/,  — .  1  =  "V  a‘+b' 

Also,  if  the  outward  normal  is  defined  as  positice  (I'igure  8),  the  relation 
local  and  global  Cartesian  CcKirdinates  is  [42] 
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Considering  (I2l)  and  using  chain  rule  of  differentation,  we  have 

>,+.t  cl'i,  a.v  d'\  dy  1,  Vi  'Vi 


(120) 

between 


(121) 


(122) 


f>7 
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^i+3  Qy.  Qs_  Qy  Qn,  ^z  1,  ',+3  1,  '',+3 


where  u^.  ^  (i=l,2,3)  denotes  the  tangential  derivatives  of  u  at  one  of  the  mid-side 

nodes  (Figure  9),  which  can  be  further  interpolated  from  the  displacements  and  their 
gradients  at  the  two  ends  of  the  corresponding  element  boundary,  i.e. 

Q  a  b  o  a  b 

u  =— -^u  — — + — ^u  +-^u.— — f-u  + — !-u  (124) 

■'1+3  21  J  41  >j  41  'j  21  41  'k  41  h 

Recalling  (lOO),  we  have 

u  =S.-T  +S,.T  (125) 

7j  N/j 

U  =S.,T  +S,  T  (I2fi) 

''k  '"k 

Substituting  (124),  (125)  and  (126)  into  (122)  and  (1230,  we  can  express  u^,  u,  at  each 
of  the  element  mid-side  nodes  in  terms  of  their  corresponding  normal  derivatives  u„,  as 
well  as  the  mixed  nodal  degrees  of  freedom  at  the  two  ends  defining  the  segment. 
That  is 
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1  i  I  i 


Following  the  same  procedure  and  considering  (103),  the  transformation  equations  for 
V  ,  Vj,  Wj.,  w^  at  each  of  the  element  mid-side  nodes  are 
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(129) 
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(132) 


Substituting  (127)-(132)  into  (119),  the  relation  between  the  surface  traction  components 
at  each  of  the  element  mid-side  nodes  and  the  corresponding  displacement  normal 
gradients  is 
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where 
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and  [Tj],  is  the  product  of  the  transformation  matrix  shown  in  (119) 
following  matrix  [TA] 
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with  the 
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(138) 

Rfiirran^inj’  (133),  the  ilisplacement-strevs  translormatons  for  the  normal  displacement 
gmdients  at  each  ol  the  mid  side  notles  in  the  l  (’(’'r  12  element  arc 
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or  symbolically, 


(140) 


where 


[II], -[f,],' 

(141) 

[1.] 

(142) 

]s}=-[T.];riv, 

(143) 

Here,  [H],  denotes  the  transformation  matrix  which  directly  related  the  normal  slope 
quantities  at  the  mid-point,  i+3,  to  the  corresponding  surface  traction  components.  [L], 
can  be  regarded  as  the  coupling  matrix  between  the  normal  gradients  and  the  mixed 
degrees  of  freedom  associated  with  the  nodal  points,  j  and  k,  at  the  two  ends  of  the 
boundary.  {S}j  is  the  local  effect  resulting  from  the  applied  uniform  loading  e„. 


5.3  Finite  Element  Formulation 

The  discretized  form  of  (7.3)  can  be  rewritten  as 
M 

=  (144) 

m=  1 

where  [K"]  is  the  element  stiffness  matrix,  {F"'}  denotes  the  nodal  force  including  the 
local  effect  due  to  the  uniform  in-plane  strain  loading  as  shown  in  (80)  for  the 
free-edge  stress  problem,  {q”}  is  the  set  of  generalized  nodal  displacement  components 
within  the  element  m,  and  M  denotes  the  total  number  of  elements. 
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In  order  to  hove  boili  displacement  and  traction  continuity  along  inlerelement  as 
well  as  interlaminar  boundaries,  assuming  (144)  is  based  on  continuous  strain  cubic- 
displacement  interpolation  as  in  IX"CI'-12  described  in  4.3.2,  the  displacement-stress 
transformation  matrices  derived  in  the  previous  section  were  imposed  on  the  generalized 
degrees  of  freedom.  Combining  (104)  and  (140),  the  generalized  displacement 
components  in  the  l.(X7r-12  element  were  transformed  to  the  mixed  type  degrees  of 
freedom  as  Inllows: 


w  here 


((i"'}-rr"'ll(ri+fi>"'l 


[(;]  0  0  0  0  0 

0  [G]  0  0  0  o 

0  0  [G]  0  0  0 

lU,  [U3  0  [ll]3  0  0 

0  [L],  [L],  0  [H],  0 

[Llj  0  [Ll  0  0  [Hlj 

(P"'}H{R}.{RMKUS}3,{S},,{S),r 


Iq"'!  =  {u  ,,v  ,.w  ,.u^^.v^  ^,w^,^.u^^,v,^.w^^.U3,v,,W3.u^^^.Vy^,w^^,u,^^^  ,w. 


u,,v,,w„u  .V  ,w  ,u  .V  ,w  .u  .\  ,w  .u  ,v  ,w  ,u  ,v  ,w  )'  (148) 

->  r  V,  /,  /,  /,  i.,,  I.J  1.4  ..4  11^  I,. 

<u'"'f=:  =  {u  ,,V  ,.W  ,,ll  .V  ,W  .T  ,T  .CT  .11  ,V  ,W  ,T  ,T  ,CT  ,U  ,,V  „W 

'  ’  '  r  I’  r  y,  V,  V,  v/|  /,  ^  :!  \y  >3’  y,  y./2  y/,’  i  <  ' 


u  ,v  ,w  ,T  ,T  ,cr  ,<r  ,tr  ,cr  ,cr  ,a  ,a  ,cr  ,cr  ,tr  }'  (149) 

''3  '3  '3  ’“'3  '"^3  '3  ""(>  '""h  "'4  "'’4  '"4  -S  ""S 

Hach  of  the  m  matrices  in  (146)  was  divided  into  two  parts  to  match  the  mixed 
degrees  of  freedom  in  (149).  In  the  finite  element  computation,  this  transformation 
was  implemented  during  the  formation  of  each  of  the  subtriangular  element  stiffness 
matrix  and  load  vector  correspmding  to  the  ). ()('!'  12  element.  Thus,  due  to 


appearance  ol  displacement  and  interlaminar  traction  components  at  the  corners  of  the 
triangle  as  well  as  traction  comjxtnents  at  the  mid-side  nodes,  a  cubic  variation  of 
displacement  and  a  quadratic  variation  of  traction  were  ensured  along  element 
boundary.  More  importantly,  both  fields  are  continuous  across  the  common  boundary 
between  two  contiguous  triangular  elements.  Substituting  (145)  into  (144),  we  have 

M 


n  =  £(l(q'"'nT'"r-t-ii>’"r)  [K"']  ([•r])q'‘''}  +  {l>’‘'})-({q'’"}'[T"’r-t-{P''’}'){F"'}) 


til  I 


.M 


tit  I 


+  l{p'"}‘[K''']{P''V{q''''r[T''f{l"''>-il>''l'n" 


(150) 


or 


a  =  £(  i{q'"’r(K'"]{q''"}-{q''"r{R'"}+C'") 


Tn»  I 

where 

[i<’"]=[T'"f  [K"’]  [T"’] 
{R"'}=[T'"f{F"'HT"'f[K '"]{?"’} 

C:"'  =  l{p'"nK"‘]{l>"'}-<P"‘r<F"')=constant 

Using  global  cixirdinates,  (151)  could  be  written  as 


n=l{qf[K]{q'Hq’}m}+C„ 


(151) 

(152) 

(153) 

(154) 

(155) 


Taking  the  variation  of  (155)  with  respect  to  iq’}  yields  the  system  generalized  nodal 
force-displacement  equilibrium  equations 

m{q'}={m  (156) 
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'I'lie  displacement  as  well  as  interlaminar  stress  comixments  were  then  obtained  by 
solving  the  resulting  set  of  linear  equations  in  the  standard  manner. 

5.4  Calculation  of  Stresses 

Solutions  of  the  finite  element  system  gives  the  three  displacement  components, 
their  tangential  gradients  along  the  element  edges,  rotation  about  the  longitudinal  axis 
and  interlaminar  stress  field  at  the  tinner  iumJcs  of  the  Q  23  element,  along  with  the 
Unindarv  traction  conifKinents  at  the  Lenter  point  o)  each  of  the  Q-2.^  element  surfaces. 
To  retrieve  the  rest  of  the  displacement  components  at  each  corner  node,  the 
transformation  matrix  used  in  (104)  was  reapplied  to  the  calculated  nodal  point 
solution.  Having  found  the  complete  displacement  gradients  and  hence  strains,  the 
inplane  stress  components  at  the  corner  nodes  of  the  Q-23  element  could  then  be 
computed  using  (115). 

The  interlaminar  or  interelement  stress  components  at  the  mid-point  of  each 
element  boundary  are  merely  the  traction  components  directly  produced  by  the  solution 
of  the  continuous  traction  Q-23  element  provided  a  rectangular  mesh  is  used  in  the 
analysis.  For  determining  the  rest  of  the  stress  components,  the  normal  displacement 
gradients  at  mid-side  ncxles  were  recovered  from  the  mxlal  point  solution  using  (140). 
Furthermore,  the  displacement  and  their  gradients  along  the  edge  at  the  mid-point  were 
interpolated  from  the  previously  computed  displacement  components  at  the  two  ends  of 
the  segment.  Having  transformed  these  displacement  gradients  from  local  to  global 
Cartesian  coordinates,  calculation  of  the  remaining  strain  and  stress  components  at  the 
mid-side  of  each  of  the  element  boundaries  is  direct. 
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5.5  Boundary  Conditions  of  a  Quadrant  oi"  the  Delamination 
Specimen 

For  the  free-edge  delamination  specimen,  because  ol'  symmetries  in  the  laminate, 
only  one  quadrant  of  an  x=constant  plane  was  considered  (Figure  2).  Along  the 
boundary,  either  displacements  or  tractions  are  specified  at  each  point. 

5.5.1  Boundary  Conditions  Along  Lines  of  Symmetry 


Symmetry  of  loading  as  well  as  geometry  alxiut  the  mid-plane  implies  that  the 
displacement  functions  satisl'y  the  following  condition 


L(y.-/.)=U(y,zl 

(157) 

\'(y,-z.J=V(y.z) 

(158) 

W(y,-z)=-W(y,z) 

(159) 

Using  chain  rule  of  differentiation, 

-U/y,-z)=U^(y,z) 

(160) 

— V^(y,-z)=V^(y,z) 

(161) 

Wy(y,-z)=-Wy(y,z) 

(162) 

Setting  z=0  in  (5.79)-(5.82),  we  have 

W(y, ())=() 

(16.^) 

W  (v,0)=U  (y,())=V  (v.())={) 

(164) 

From  (5.84),  >„(y,0)=yj./y,0)=0,  and  consequently,  for  the  layered  orthotropic 

material. 

(165) 

Also 

(w^-v^Xy.o)  =  0 
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Invariance  under  a  rotation  ol'  180  degree  alrout  the  z-axis  through  the  center  of  the 
specimen  implies 


L'(-y,z)=-lj(y,z)  (166) 

V(-y^)=-V(y;z)  (167) 

W(-y^)=W(y,z)  (168) 

(166)  and  (167)  lead  immediately  to 

L((),/.)=\'((),z)=()  (169) 


lor  all  /..  and  conse(|uentl v 

r^((),/)=\ ^(().z)-()  (170) 

llv  chain  rule  of  dil  ferentiation,  (166)  yields 

-Wy(-y,z)=W^(y,z)  (l7l) 

Hence,  for  y=0 

W^(0,z)=0  (172) 

(170)  and  (172)  imply  •y„(0,z)=Vyj(0,z)=0.  Hence,  for  the  layered  orthotropic  material, 
T^/0,z)=Ty/0,z)=0  (173) 

Also 

(\V  -V  /0./,):::0 

(,'omhining  (163).  (164),  (165),  (169).  (172)  and  (173),  along  with  L(0.0)=0  in  order 
to  prevent  rigid-body  displacement  of  the  laminate,  the  continuous  traction  finite 

element  model  for  a  quadrant  laminate  under  consideration  should  satisfy  the  following 


conditions  along  lines  of  symmetry 

U((),0)=U(0,z)=V(0,z)=W(y,0)=W^(y.0)=\V^(0,z)=0 


(174) 


5.5.2  Traction-Free  Boundary  Conditions 

The  traction  tree  Ix'undary  condiiions  associated  witli  one-quadrant  of  the  laminated 
specimen  are 

T^^(y,H)=Ty^(y,H)=o-^(y,H)=0  (176) 

at  the  lop  surface,  along  with 

cr  (I5,z)=t  (n./.)=T  (n,z)=0  (177) 

\  \\  V/ 


at  the  lateral  I  rce  edge.  Here  211  and  211  denote  the  total  width  and  thickness 
1  esjX'et i\ el V  ol  the  laminated  s]xcimen.  '  sing  the  continuous  trtiction  f,)  23  model, 
irtiction  free  bound. ir\  conditions  slu'W  n  in  (176)  c.m  be  identically  satisfied  lor  nod.d 
points  t>n  the  top  surface.  However,  due  U'  the  lack  of  inplane  stress  comjxments  as 
ntxlal  degrees  ol'  freedom,  only  the  last  condition  shown  in  (177)  can  be  s}-)ecified  at 
those  element  corner  nodes  along  the  lateral  free-edge.  To  completely  satisfy  the 
traction-free  condition  along  the  lateral  free-edge,  the  following  device  was  developed. 

In  order  to  enforce  the  remaining  two  inplane  stress-free  conditions  in  (177),  it  is 
necessary  to  express  these  in  terms  of  nodal  point  degrees  of  freedom.  This  results  in 
a  linear  relationship  lietween  the  degree.s  of  freedom  at  nodal  points  on  the  free  edge. 

The  stress  free  condition  can  be  written  explicitly  as 


a  =<)=r,,e  +(T  \v  +('„  u 

V  le  o  -1  12  V  !(•  V 


T  =()=(',  e  w  V  +(',  ,u 

\  1(,  ll.  /  1(,  V  M.  V 

or  symbolically 


(17S) 
( 1  79 ) 


ic  c  c  r 

'  12  23  22  26 

r  r  c  c 

16  36  ^26  66 


(180) 
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\\'here  are  componenis  of  tlie  slilf'ness  matrix  del'ined  by  (116).  Solving  above 
equation,  the  inplane  strain  components  u,  v  can  be  expres.sed  in  terms  of 
interlaminar  normal  strain  by 

(181) 

i^y)  h. 

where 

I  =1.(6:  r  -r  r  j 

'll  ^  .'1.  II,  l.l.  1-'^ 


V  1 
y 

*11 

*12 

e 

c 

u 

yl 

*21 

*22 

W 

1  Z 

and 


(182) 


(183) 


Substituting  (102)  into  (I8I)  for  w^,  we  can  further  relate  u^.  and  v^,  to  the 
interlaminar  normal  stress  component  tr^  through  the  following  linear  relationships 

v,=q|Cr^+q,e  (185) 

w  here 


I  I *22®2~^  12*21^2^*22®! "*^^21^ 

qi=jti,A -X)] 


^2  =  -^(',2':.,**-*i.':2l*,  +  'i2**,+*,i* 


(186) 

(187) 

(188) 

(186) 
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and 


/3=1-1. 


(190) 


Again,  li,,  B^,  B^  and  X  occurring  in  (186)  through  (190)  have  been  defined  by 
(96)-(99).  Thus,  the  variables  associated  with  the  nodal  points  of  the  free-edge  in  the 
lateral  continuous  traction  Q-23  element  are  no  longer  independent  but  related  through 
(184)  and  (185).  Incorporating  these  linear  relationships  in  the  displacement-stress 
transformations  slunc  n  in  (104)  and  (14(0  fvir  elemenis  on  the  lateral  boundary  implies 


satisfaction  ol'  tlie  traction  free  houndarv  i(>nditions  (177).  The  transformation  becomes 


u 

1  0  0  (1  0 

0 

0 

0 

0 

11 

0 

\ 

0  1  0  0  0 

0 

0 

0 

0 

\ 

0 

w 

0  0  1  0  0 

0 

0 

0 

0 

w 

0 

u 

>• 

0  0  0  0  0 

0 

0 

0 

Pj 

u 

y 

P2^> 

V' 

V 

= 

0  0  0  0  0 

0 

0 

0 

‘I. 

V 

y 

+ 

w 

y 

0  0  0  0  0 

1 

0 

0 

0 

w 

y 

0 

“z 

0  0  0  0  0 

0 

§4.5 

0 

7 

xz 

0 

V 

z 

0  0  0  0  0 

-1 

§45 

§44 

0 

T 

yz 

0 

W 

z 

0  0  0  0  0 

0 

0 

0 

rj 

cr 

z 

2  o 

where 

r2  =  B,+Bjq^+B,p, 


(192) 

(19.1) 
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SECTION  VI 


ANALYSIS  OF  FREE-EDGE  DELAMINATION  IN 
LAMINATE  COMPOSITE  SPEOMENS 

6.1  Introduction 

Continuous  traction  finite  element  formulation  developed  in  the  previous  section 
was  applied  to  obtain  an  approximation  to  displacement  as  well  as  stress  fields  in  a 
free-edge  delamination  specimen.  The  analysis  consisted  of  two  parts.  The  first 
consisted  of  solving  four-ply  symmetric  laminates.  The  purpose  was  to  examine  the 
credibility  of  the  continuous  traction  finite  element  model  by  comparing  the  numerical 
solutions  with  those  from  Pagano's  analysis  based  on  a  generalization  of  Reissner's 
theory.  At  the  same  time,  the  deficiency  in  using  the  continuous  strain  free-edge 
stress  model  described  in  Section  IV  was  examined.  The  second  part  consisted  of 
studying  of  edge  delamination  tendency  in  two  classes  of  multi-ply  laminates  subjected 
to  longitudinal  loading.  For  laminate  specimens  with  stacking  sequence  of 
[(0/— 0)^/9O„  j],,  displacement  and  stress  fields  calculated  from  the  continuous  traction 
Q-23  element  were  compared  with  those  obtained  using  an  overlay  procedure  with 
constant  strain  elements  [45]  and  with  experimental  observation  [46].  For  the 
[(©/— 0)^90],  laminate  specimens,  the  continuous  traction  finite  element  code  was  used 
in  conjunction  with  some  well  known  anisotropic  failure  criteria  [47-50]  to  predict  the 
onset  of  transverse  cracking  and  the  onset  of  edge  delamination  in  various  laminate 
specimens  observed  in  experimental  data  [51].  The  purpose  was  to  evaluate  the 


suitability  ol  these  I'ailure  (.riteria  lor  application  to  the  free-edge  delamination  problem. 
Ikcause  of  symmetries  in  the  laminates,  only  one  ipiadrant  (Ihgure  2)  was  considered 
in  each  case. 

6^  Four-Ply  Laminates 

In  this  .section,  analy.sis  of  two  long  .symmetric  laminate  strips  made  of 
graphite-ep»i\y  materials,  with  I'ilrer  orientations  of  [45-' 45]  and  [()/f>()],  under  uniform 
inplane  strain  in  the  longitudinal  ilirection  is  descrilxfii.  I'he  relation  Iretw'een  laminate 
width  and  thickness  was  2b^l6h  lollowing  [8].  In  the  analysis,  each  ply  was 
idealized  as  a  homogeneous,  elastic  orthotropic  material.  Tor  comparison  purpose,  the 
material  properties  a.ssumed  here  following  Pagano's  work  [8] 

E,,  =20X10^1 

E22=Ejj=2.1Xl0^psi 

G  j  2=G  I  j=G2j=0.85x  lO^psi 

‘'i2=*'i3=>^23=0.21 

The  subscripts  1,  2  and  3  correspond  to  the  longitudinal  transverse  and  thickness 
directions  respectively.  A  144-element  model  as  shown  in  Figure  (lOa)  was  used  to 
discretize  a  typical  x=const;mt  plane.  Numerical  results  based  on  the  continuous 
traction  Q-23  and  continuous  strain  (,1-23  elements  were  compared  with  Pagano's  [S] 
analytical  solution. 

Figures  (11)  through  (23)  illustrate  comparisons  for  both  stress  and  displacement 
fields  at  specific  locations  for  the  angle-ply  and  cross-ply  laminates  using  different 
solution  techniques.  The  value  of  N  in  these  figures  corresponds  to  the  number  of 
sublayers  used  in  Pagano's  theory.  Thus,  N=6  indicates  that  each  physical  layer  of 
thickness  h  was  modeled  by  three  sublayers  each  of  thickness  h/3,  wdrile  N=2  denotes 
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that  each  physical  layer  is  treated  as  a  unit  [S].  Alsn.  the  calculated  displacements  and 
stresses  were  normalized  by  the  applied  uniform  strain  loading  e„  which  has  been 
taken  as  unit  in  the  present  analysis. 

6.2.1  Angle-Ply  Laminate 

Figures  (ll)  and  (12)  show  the  distribution  of  cr^  and  along  the  width  of  the 
laminate  at  the  center  line  of  the  top  (45  degree)  layer.  The  results  obtained  using 
the  continuous  traction  element  agreed  (juite  well  with  those  of  Ptigano's  \-h 

solutions  across  the  entire  width  ol  the  laminate. 

A  comparison  ol'  the  shear  stress  (t^^)  distribution  along  the  interface  ol  the 
45.-45  layers  (ligure  1.^),  indicated  that  the  continuous  traction  Q-23  solution  had 
sharp  ri.se  toward  the  free-edge  similar  to  Pagano's  solution  with  N=6.  Satisfactory 
agreement  was  observed  between  these  two  solutions  for  stress  across  the  width  except 
at  the  free-edge  boundary  where  continuous  traction  Q-2.3  somewhat  underestimated  the 
singular  stress.  Figure  (14)  shows  the  through-thickness  stress  distribution  of 
calculated  from  both  continuous  strain  and  continuous  traction  0-23  elements  at  the 
free-edge  of  the  laminate.  Very  close  agreement  was  generally  observed  between  these 
two  solutions  throughout  the  thicknes.s.  Also,  at  the  interface  of  the  45/-45  layers, 
continuity  of  the  interlaminar  shear  stress  w'as  ensured  for  both  cases.  This  is  because 
of  the  rotational  symmetry  about  z-axis  in  the  particular  angle-ply  laminate  considered. 
The  singular  behavior  of  which  is  highly  localized  at  the  interface  between  45/-45 
layers,  is  noticeable.  The  distribution  of  along  the  interface  between  45°  and  -45“ 
plies,  which  was  not  indicated  in  [8],  will  be  discussed  in  the  next  section. 

For  the  axial  displacement  distribution  across  the  w'idth  of  the  top  surface, 
continuous  traction  Q-23  results  compared  w'ell  with  Pagano's  N=6  solution  (Figure  15). 
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Figure  15: 


Axial  Displacement  Across  'I'op  Suriace  of  [45/-45|  laminate 


I'igure  (16)  slio^^s  the  thniugh-thickness  (.lisirihutioti  ol'  ;ixia)  displacement  based  on 
continuous  strain  as  well  as  continuous  traction  (J-23  elements  at  the  free-edge  of'  the 
laminate.  Again,  these  two  .solutions  matched  well  throughout. 

6.2.2  Cross-Ply  Laminate 

Distribution  of  cr,  along  the  width  on  the  central  plane  of  the  [0/90l  laminate, 
show  n  in  l  igure  (17),  indicates  a  sharp  rise  near  the  f ret -edge  boundary.  .Solutions 
obtaineil  from  the  contiiuu'us  traclum  Q  2^  elemcni  nearly  coincided  with  Pagano's  \  (> 
solution  (wer  the  entire  width  of  the  laminate. 

figure  (18)  sliows  the  variation  of  ct,  along  the  interface  lietw'een  the  O'  and  9o 
plies.  Due  to  the  presence  of  the  discontinuity  in  elastic  properties,  a  singular  stress 
behavior  w'ould  be  expected  at  the  free  edge.  On  the  boundary,  re.sult  from  the 
continuous  traction  Q-23  element  had  a  steeper  gradient  than  that  of  Pagano's  theory. 
Apparently,  one  possible  reason  for  this  discrepancy  is  that,  in  Pagano's  analysis,  each 
physical  layer  could  be  modeled  by  at  most  three  sub-layers.  However,  in  the  finite 
element  analysis,  the  thickness  was  divided  into  eight  elements.  If  a  coarser 
discretization  were  to  be  used,  say  4x18,  cr^  calculated  from  the  continuou  traction 
Q  2^  element  would  possibly  agree  quite  well  with  Pagano's  solution  at  the  free  edge 
interface.  figure  (19)  shows  the  influence  of  through  the  thickness  refinement  of 
mesh  on  cr,.  figure  (20)  shows  through  ihe-thickness  distribution  of  ct,  at  the 
free-edge  of  the  laminate.  In  the  vicinity  ol  the  interface,  the  continuous  traction 
Q-23  element,  enforcing  continuity  of  cr,,  at  the  interface  betw'een  differently  oriented 
layers,  gave  a  stress  distribution  quite  different  from  that  given  by  the  continuous 
strain  Q-23  analysis.  Away  from  the  interface  the  two  sets  of  results  were  close. 
Also,  the  interlaminar  stress  cr,  was  observed  to  ha\e  a  maximum  value  in  the  interior 
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DISTANCE  FROM  CENTER  LINE  Y/B 


Figure  17:  Distribution  of  Z-stress  Along  Central  Plane  of  [0/90]j Laminate 
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Figure  19:  Effect  of  Mesli  Refinement  on  tlie  Z-stress  Distribution  Along 

0/90  Interface 
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Figure  20:  Through-Thickness  Distribution  of  Z-stresss  at  the  Free-Edge  of 

[0/90;^  Laminate 
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of  the  ’^O-deg  layer  closer  to  the  interface  with  the  top  layer.  However,  in  both  cases, 
the  solutions  displayed  oscillatory  patterns  near  the  interlace.  This  could  be  due  to  the 
finite  element  mesh  used  being  not  fine  enough  to  approximate  the  steeply  varying 
stresses  associated  with  abrupt  change  of  material  properties. 

Values  of  t^,,  along  the  interface  between  the  [0/90]^  layers,  calculated  from  the 
continuous  traction  (^-2,1  element  (1-igure  21),  showed  satisfactory  agreement  with  those 
calculated  by  Pagano.  This  is  Irecause  the  continuous  traction  Q-23  element  exactly 

satisfies  the  traction-free  boundtiry  condition  similar  to  Pagano's  theory,  llowexer.  an 

oscillatory  error  was  obserxed  near  the  free  edge.  Apparently,  further  mesh  refinement 
along  the  y-direction  is  required  near  the  free-edge  in  order  to  approximate  the  singular 
stress  behavior.  IMgure  (22)  displays  through-the-thickness  stress  distribution  ol’  r,., 
calculated  from  both  continuous  strain  and  continuous  traction  0-23  elements  at  the 

free-edge  of  the  laminate.  Apparently,  satisfaction  of  the  traction-free  boundary 
condition  associated  with  the  continuous  traction  0‘23  element  represents  an 

improvement  over  the  continuous  strain  0-23  element. 

Comparative  results  for  the  variation  of  transverse  displacement  along  the  top 
surface  of  the  [0/90]..  laminate  are  showm  in  Figure  (23).  Excellent  agreement  was 
observed  between  results  using  the  continuous  traction  0'23  element  and  Pagano's  N=2 
solution. 
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6.2.3  Effect  of  Traction-Free  Edge  on  the  Solutions 

In  order  to  investigate  the  effect  ol'  requiring  tlie  satisfaction  of  a  traction-free 
boundary  condition  on  the  finite  element  solutions,  the  continuous  traction  Q-23 
element  was  employed  with  only  the  requirement  that  along  the  lateral 

free-edge  of  the  four-ply  laminate  specimens.  In  other  words,  the  displacement 
constraint  conditions  developed  in  (188)  and  (189)  used  to  specify  the  in-plane 
stress-free  Iroundary  conditions  were  not  impi'sed  in  this  model,  for  convenience  in  the 
following  comparisons,  this  is  designated  as  continuous  traction  (partial). 

I'igure  (24)  shows  the  distribution  of  along  tlie  interface  of  the  45/-45  layers. 
Solutions  calculated  from  the  continuous  traction  (partial)  hud  a  steeper  gradient  in  the 
vicinity  of  the  free-edge  than  that  of  previous  continuous  traction  0"23  element.  A 

similar  observation  is  made  for  the  variation  of  cr^  at  the  interface  of  [0/90]^  laminate 

(Figure  25).  Thus,  it  is  concluded  that  the  nonimposition  of  the  conditions 
cry=0  and  Tyj=0  would  overestimate  the  magnitudes  of  the  interlaminar  stresses  on  the 
interface  near  the  free-edge.  However,  the  nonsatisfaction  of  these  two  traction-free 
boundary  conditions  had  no  significant  effect  on  the  displacement  field  in  the  laminate 
specimen.s.  Figure  (26)  shows  the  solutions,  for  axial  displacement  distribution  across 

the  width  of  the  top  surface  in  the  (45'-45].  laminate,  obtained  from  the  continuous 
traction  (partial)  and  from  the  continuous  traction  Q-23  element  which  satisfies  all 
traction  free  conditions  at  the  free-edge. 
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Figure  24: 


Distribution  of  XZ-stress  Along  i45  Interface— Effect  of 

Traction-Free-Edge 


Z-STRESS/^.X  10"'’  (PSn 
.10  0,00  0,  10  0,20  0,30  0,40 


0  COMPLETE  TRflCTION-FREE 
5^  PARTIAL  TRACTION-FREE 


ri - i - 1 - 1 - 1 - 1 

0.00  0.20  0.40  0.60  0.80  1.00 

DISTANCE  FROM  CENTER  LINE  Y/B 


Figure  25:  Distribution  of  Z-stress  Along  0/90  Interface— Effect  of  Traction- 

Free-Edge 
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6^.4  Effect  of  Mesh  Refinement 

The  analysis  of  the  four-ply  laminate  specimens  was  originally  carried  out  by 
using  the  continuous  traction  Q-23  element  w'ith  a  uniform  80-element  model  as  shown 
in  Figure  (27).  However,  for  reliability  of  results,  the  finite  element  mesh  must  be 
refined  in  regions  of  steeply  varying  stresses.  This  results  in  the  present  analysis 
using  the  144-element  model  w'hich  indeed  was  obtained  by  dividing  the  two  elements 
closest  to  the  free-edge  in  the  80-element  model  into  10  elements  along  the  y-direction. 
To  study  the  effect  of  mesh  refinement  on  the  continuous  traction  finite  element 
solutions,  comparisons  were  made  for  the  stress  distributions  in  the  four-ply  laminate 
specimens  between  the  uniform  80-element  and  the  kx;ally  refined  144-element  mixlels. 

Figure  (28)  shows  the  distribution  of  at  the  interface  of  [45/-451  laminate 
ba.sed  on  the  144-element  model.  A  steeper  gradient  of  was  observed  on  the 
boundary  as  compared  with  the  result  using  80-element  model.  A  comparison  of 
distribution  along  the  interface  of  the  45/-45  layers,  indicates  that  the  144-element 
model  had  a  compressive  finite  maximum  value  at  the  free-edge  rither  than  a  tensile 
quantity  from  the  80-element  model  (Figure  29).  This  indeed  has  demonstrated  the 
inappropriate  sign  of  o’,  shown  in  Figure  (1)  based  on  the  perturbation  technique  as 
well  as  finite  difference  method.  For  variation  of  o’,  along  the  interface  of  [0/90]^ 
laminate.  Figure  (30)  indicates  that  a  singular  stress  behavior  was  properly  reproduced 
by  the  144-element  model  and  did  not  show  well  in  the  results  based  on  the 
80-element  mesh. 

Use  of  144-element  model  might  still  be  insufficient  to  approximate  the  singular 
stress  behavior.  One  example  is  the  7^  distribution,  which  had  an  oscillatory  pattern 
near  the  free-edge  along  the  interface  of  [0/90],  laminate,  as  mentioned  before.  To 
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overcome  this,  a  more  rel'ined  mesh  (2()8-element  model)  obtained  by  further  dividing 
the  two  elements  closest  to  the  free-edge  in  the  144-element  model  into  10  elements, 
was  used  in  the  analysis.  Along  with  the  results  calculated  from  80-  and  144-element 
models  respectively.  Figure  (31)  indicated  the  improvement  of  the  distribution  along 
0/90  interface  over  the  boundary  layer  region  as  more  refined  elements  were  used  near 
the  free-edge.  However,  regardless  of  mesh  patterns  used,  there  was  an  oscillatory 
error  in  r,^  in  the  two  clenienis  next  to  the  free-edge. 

I'igure  (32)  slioxxs  through- ihe-thickness  distributions  of  cr^  at  the  free-edge  of 
[0'90],  laminate  based  on  the  144-element  mtxiel  but  with  finer  mesh  near  the 
interface  (Figure  10b)  and  its  refinement  (28S-element  model)  in  the  thickne.ss  direction. 
It  is  oltserved  that  the  oscillatory  error  near  tlie  interface  w'as  reduced  by  using  the 
refined  144-element  model  and  was  nearly  disappeared  under  more  refinement  over  the 
laminate  thickness.  Jvleanwhile,  the  maximum  value  of  cr,  within  the  90-degree  layer 
was  moving  clo-ser  to  the  0/90  interface. 
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Figure  31:  Distribution  of  YZ-stress  Along  0/90  Interface— Effect  of  Mesh 
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6.3  Free-Edge  Delamination  in  Multi-Ply  Laminate  Specimens 

Analysis  of  the  four-ply  laminate  specimens  described  in  the  previous  section 
demonstrated  some  validity  of  using  the  proposed  finite  element  procedures  in  solving 
free-edge  effect  problems.  Both  continuous  strain  and  continuous  traction  Q-23  elements 
had  similar  prediction  on  the  displacements  and  inplane  stress  distributions,  which  also 
compared  well  with  Pagano's  analytical  solutions.  However,  discrepancy  between 
continuous  strain  and  continuous  traction  finite  element  models  was  apparent  for  tiie 
interlaminar  stresses  near  laminate  interfaces  l:)etween  difl'erently  oriented  layers  or  near 
the  traction-l'ree  Ixiundary.  Due  to  the  fact  that  stress  continuity  across  interlaminar 
boundary  as  well  as  traction-free  boundary  condition  are  exactly  stitisfied  in  the 
continuous  traction  Q-23  element,  solutions  obtained  from  this  approach  were  expected 
to  be  more  reliable  than  those  from  the  continuous  strain  Q-23  element.  In  this 
section,  application  of  the  continuous  traction  Q-23  element  to  investigate  the  free-edge 
effect  as  well  as  initiation  of  edge  delamination  in  the  multi-ply  laminates  is  described. 

6.3.1  Analysis  of  [(6/-6)ni/90y2]sLan^i*^^tes 

Four  types  of  laminates  with  predetermined  fiber  orientations  [46]  were  used  in 


present 

investigation.  These 

are 

Type 

Stacking  Sequence 

Width 

Ply  thickness 

Plies 

A 

[(49.8/-49.8)^/90]^ 

1.0  in 

0.00506  in 

22 

B 

[(30.87-30.8)5/90]^ 

1.0  in 

0.00508  in 

22 

C 

[(25.57-25.5)5/90]^ 

1.0  in 

0.00505  in 

22 

D 

[(47.97-47.9), o/90]^ 

1.0  in 

0.00499  in 

42 

The  material  used  in  the  study  was  AS4/3501-6,  graphite-epoxy,  and  the  elastic 
constants  were  [46] 
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i:,,  =  l9.26Xl0‘’psi 
l;32=1-32xl()‘‘psi 

G,^=0.83xl()'‘psi 

2=0.35 

All  the  specimens  have  been  investigated  both  analytically  and  experimentally  at 
the  AFWA1./AFFUL  [46].  The  I3elamination  Moment  Coefficients  (DMC)  were  derived 
and  used  to  evaluate  quantitatively  the  delamination  tendency  of  the  laminates.  Alsu, 
generalized  constant  strain  element  was  applied  to  analyze  half  of  the  width  of  the 
laminate  specimens.  In  the  experimental  as|')ect,  \arious  techniques  including  I’ransverse 
Strain  Gages,  (Tracked  Silver  Ink  Instrumentation  and  Acoustic  Emission  Instrumentation, 
etc.  were  used  to  determine  the  onset  of  delamination  and  to  validate  the  analytical 
results. 

6.3. 1.1  Numerical  Evaluation 

A  154-element  model  shown  in  Figure  33(a)  was  used  to  discretize  one  quadrant 
of  a  typical  x=con.stant  plane  in  the  laminate  specimen  Types  A,  B.  C,  and  a 
294-element  model  (Figure  33(b))  w'as  used  for  speicmen  Type  D.  Each  ply  was 
modeled  by  a  single  element  through  its  thickness.  Interlaminar  stress  field  within 
various  laminate  specimens  for  an  applied  longitudinal  average  stress  of  KXl  psi  were 
computed. 

Comparisons  of  distribution  along  the  mid-plane  of  various  multi-ply  laminates, 
(Figures  34-37)  indicate  that  the  continuous  traction  Q-23  solutions  had  sharp  rise 
toward  the  free-edge  similar  to  constant  strain  element  solutions  [46].  Satisfactory 
agreement  was  generally  observed  between  these  two  solutions  for  stresses  in  the 
vicinity  of  the  free-edge  except  on  the  boundary  where  the  stress  calculated  using  the 
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Distribution  of  Z-stress  at  Mid-plane  for  Width  of  10  Plies  of 
Specimen  C 


Q-23  eJement  was  distinctl.v  less  than  that  Irom  the  constant  strain  solution. 
However,  the  cr,  values  calculated  from  the  constant  strain  element  were  extrapolated 
from  the  interlaminar  stresses  at  z=0  obtained  by  lagrangian  interpolation  of  the  tr^ 

values  at  the  element  centroids.  This  is  unlike  the  continuous  traction  Q-23  solution 

where  cr^  at  the  free-edge  was  directly  calculated  as  nodal  degree  of  freedom  and 

would  be  expected  to  be  more  reliable. 

I'iguves  show  the  through-the  thickness  stress  distributions  of  (t^  and  r^. 

calculated  from  continuous  traction  Q  2.3  element  at  the  free-edge  of  various  laminate 
s]-)ecimens.  It  is  observed  that  for  the  same  applied  axial  stress,  specimen  D  had  the 

largest  value  ol'  normal  stress  cr^  at  the  free-edge,  followed  by  specimens  A,  B  and  ('. 
Figure  (42)  illustrates  this.  The  slope  discontinuity  of  tr,  at  the  interfaces  of  the 

free-edge  shown  in  Figures  (38)-(4l)  was  possibly  attributed  to  the  material  as  well  as 
geometrical  discontinuities  in  that  region.  Figure  (43)  indicates  that  much  smoother  or^ 
distributions  through  the  laminate  thickness  were  recovered  within  the  angle-ply 
laminae  at  a  small  distance  from  the  free-edge.  In  fact,  with  further  refinement  along 
the  free-edge,  the  solution  of  o’,  on  the  free-edge  is  even  better.  Figures  (44)-(45) 

show  the  solution  for  cr,  at  the  free-edge  and  at  y=().495  for  specimen  A  with  each 
edge  element  being  refined  into  four  elements  along  y  direction.  Also,  figure  (46) 
shows  the  functional  dependence  of  longitudinal  stress  a,  as  well  as  interlaminar 

normal  stress  O’,  on  the  fiber  orientation,  respectively,  for  the  [(0/— 6)5/90]^  laminate 

under  the  same  applied  loading.  The  ordinates  of  these  curves  are  the  respective 
values  of  o,  and  o,  at  the  intersection  of  the  mid-plane  and  traction-free  edge  of  the 
laminate.  Results  obtained  from  the  continuous  traction  Q-23  element  indicated  that 
lx)th  o^  and  o,  attained  their  maximum  values  approximately  in  the  fiber  orientation 
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0  =  30°. 

The  shear  stress  distributions  were  similar  for  these  specimens,  and  their 
magnitudes  are  relatively  smaller  than  the  maximum  normal  stress  tr^.  However,  the 

existence  of  evidently  reflected  a  defect  of  the  numerical  model  adopted  in  [46]  in 

which  the  delamination  specimen  was  treated  as  an  axially  symmetric  problem  in 
which  was  inherently  assumed  to  be  zero  throughout  the  laminate  thickness.  It 

was  noted  that  c,  distribution  had  a  sloj^e  discontinuity  at  the  mid-plane  surface 
within  the  9()-degree  layer.  I'his  does  not  appear  to  be  reasonable  for  the  present 
symmetric  laminate  specimens.  Presuming  that  this  was  as.sociated  with  the  use  of  a 
single  element  through  the  thickness  of  90-degree  layer  Ireing  insufficient  to 

accommodate  the  mismatch  at  the  interface  between  angle-ply  and  cross-ply  laminae,  a 
study  was  carried  out  refining  the  mesh  near  the  midplane.  Figure  (47)  shows  the 
dramatical  improvement  of  cr,  distribution  near  the  mid-plane  surface  of  specimen  A 
as  increasing  number  of  elements  was  used  in  the  discretization  of  90sJegree  layer.  At 
the  same  time,  the  maximum  value  of  <r^  in  the  interior  of  the  transverse  layer  was 
observed  to  move  closer  to  the  interface  with  the  angle-ply  layer.  Again,  if  both  the 
90-degree  layer  and  the  free-edge  elements  were  refined,  the  improvement  of  cr,  was 
not  only  on  the  90-degree  layer  but  also  on  the  entire  laminate.  Figures  (48K49) 
show  this  improvement. 

Comparison  of  the  interlaminar  shear  stress  Ty^  at  the  center  line  of  90-degree 
layer  along  the  width  of  various  laminate  specimens  are  shown  in  Figures  (50)-(53). 
Solutions  from  both  methcxls  indicate  that  approached  finite  maximum  values  in  the 
vicinity  of  free-edge  yet  the  traction-free  boundary  condition  could  only  be  satisfied  by 
using  the  continuous  traction  Q-23  element.  The  maximum  values  of  from  the 
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Figure  41:  Through-Thickness  Stress  Distributions  at  the  Free-Edge  of  Speci¬ 

men  D  Based  on  the  Continuous  Traction  0-23  Element 
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Figure  53;  Distribution  of  YZ-stress  at  the  Center  of  90-degree  Layer  for 
Width  of  8  Plies  of  Specimen  D 
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Figure  56: 


Through-Thickness  Distribution  of  YZ-stress  at  the  Center  of  the 
Second  Flement  From  the  Free-l-,dge  for  Specimen  C 
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two  models  were  compar;ih!e.  The  constant  strain  triangular  element  approximation 
departed  significantly  from  the  Q-23  solution  in  the  vicinity  of  the  free-edge.  This 
could  largely  be  due  to  the  nonsatisfaction  of  the  traction-free  boundary  condition 
Figures  (54)-(57)  illustrate  through-the-thickness  stress  distribution  of  along  center 
line  of  the  second  element  from  the  free-edge  for  various  laminate  specimens.  1'he 
reason  to  choose  this  site  for  comparison  was  because  of  the  singular  stress  behavior 
being  displayed  in  the  finite  element  discretization  for  both  methods  (Figures  5()-5,f). 
A  close  agreement  was  generally  observed  between  these  two  solutions.  The  continuous 
traction  (J-2.^  analysis  show  that  the  maximum  occurred  at  the  interface  Iretween 
the  negative  angle-ply  and  the  9()-degree  layers  for  all  the  specimens.  The  constant 
strain  element  does  not  have  the  capability  to  predict  this.  It  is  also  noted  that  for 
the  same  applied  axial  stres.s,  specimen  D  has  the  largest  value  of  followed  by  A, 
B  and  C,  similar  to  the  observation  for  cr^. 

Table  (2)  shows  the  values  of  interlaminar  normal  stress  cr^  at  the  interface  of 
the  free-edge  for  laminate  specimens  A,  B,  C  and  D  based  on  constant  strain  element 
and  continuous  traction  Q-23  element,  respectively.  The  ratios  of  the  normal  stresses 
CTj  to  the  relatively  minimal  value  among  them  are  2.11  :  1.32  :  1.0  :  2.47  for 
constant  strain  element,  and  2.12  :  1.33  :  1.0  :  2.39  for  continuous  traction  Q-23 
element.  Thu.s,  the  normal  stress  ratios  for  specimens  A,  B,  C  and  D,  calculated  from 
these  two  finite  element  schemes  are  comparable  with  each  other. 

Figure  (58)  shows  exaggerated  views  of  the  displacement  fields  based  on  the 
continuous  traction  Q-23  element  in  specimens  A,  B  and  C.  Figure  (59)  shows  the 
distortion  of  Specimen  D.  The  maximum  displacement  in  the  y-direction  calculated 
from  both  constant  strain  element  and  continuous  traction  Q-23  element  shown  in 
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'J'able  2:  ('omparisons  ol  Delamination  'I'endency  for  Various  Speimens 

(  unit  .  psi) 


Classification 

A 

B 

C 

D 

cr^  from  Constant 
Strain  element 

Value 

107 . 74 

67.65 

51.12 

126 . 22 

Ratio 

2.11 

1.32 

1.00 

2.47 

1 

cr  from  Continuous 
/ 

Stress  Element 

Value 

82  07 

51.54 

38.68 

- 1 

92.43 

Ratio 

2  12 

1 .33 

1.00 

2.39 

a 

Note  ;  Ratio= - — 

cr 

^inin 
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(a)  Undeformed  Model  of  D  laminate 


(b)  Deformed  Model  of  D  Laminate 


Figure  59; 


Finite  Element  Models  D  Lamiiiate  Specimen 


'I'ahle  (3),  were  found  to  Ire  jn  reasonable  agreement  for  all  these  specimens.  Also, 
from  these  figures,  it  is  observed  that  under  the  siime  applied  axial  stress,  specimens  D 
had  the  largest  distortion  near  the  free-edge,  followed  by  laminates  A,  B  and  C  in 
descending  order.  In  other  words,  specimen  D  had  the  greatest  tendency  for  edge 
delamination  among  the  four  specimens,  and  the  delamination  would  set  in  at  the 
lowest  applied  axial  stress.  This  again  confirms  the  prediction  based  on  the 
interlaminar  norma!  stress  cr  which  had  the  largest  values  for  specimen  D  as  shown 
in  Table  (2).  Ihised  u|X)n  alxive  analysis,  we  conclude  that  both  the  distortions  and 
the  values  of  normal  stress  a,  near  the  Iree-edge  of  the  specimens  calculated  I'rom  the 
continuc'us  traction  Q-23  element  are  consistent  with  those  of  the  constant  strain 
element  [46]  which  is  much  more  economical  to  implement. 
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Table  3;  ('omparisons  ol  Maximum  'I'ransverse  Deformation  for  Various  Specimens 


(  10  in) 


Classification  Spe  cimen 


Constant 
Strain  Element 


Continuous 
Traction  Element 


Value 

0 . 1309 

0.1032 

0.0829 

0.1935 

Ratio 

1  5789 

1.2448 

1.00 

2.3337 

Value 

0,1171 

0  0888 

0.0704 

0.1638 

Ratio 

1  66-.  1 

1.2619 

1  00 

2.3274 

Mote  :  Ratio=- 
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6.3.1^  Analytical-Experimental  Correlation 


An  experimenta]  study  was  conducted  at  Al'WAL/AFFDL  [46]  to  validate  the 
analytical  results.  Part  of  the  test  results  are  summarized  in  Table  (4),  which  gives 
the  specimen  type,  strains,  and  stresses  for  the  initiation  of  edge  delamination 

determined  by  transverse  gages,  silver  ink,  acoustic  emission  and  visual  observation. 

The  table  also  contains  average  axial  stress  I'or  initiation  of  edge  delamination 

calculated  from  the  continuous  traction  Q-23  element.  ’I'he  average  initial  delamination 
stresses  for  laminates  A.  II.  and  C  were  found  to  be  19.2.  23.3  and  26.4  (ksi). 
respectively  [45].  The  corresponding  finite  element  solutions  were  15.883,  25.817  and 
29.935  (ksi).  To  calculate  these  values,  the  nodal  axial  stresses  tvithin  each  LC(rr-12 
triangular  element  had  to  be  recovered  from  the  stresses  calculated  along  the  boundary 
of  (>23  quadrilateral.  The  axial  loading  applied  on  each  triangular  element  could  then 
be  computed  by  integrating  the  nodal  stress  values  over  each  triangular  area.  Having 
assembled  the  element  axial  loading  over  the  whole  system,  the  average  axial  stress 
was  obtained  by  dividing  the  total  axial  loading  by  the  total  cross-sectional  area. 


Table  4:  lixperimental  Results  Ibr  Various  luiminate  Specimens 


Laminate  Specimen 

A 

B 

1 - 1 

Axial  Strain 

10'*’ 

6099 

3300 

2900 

Initial  Matrix  Cracking  Stress 
(  ksi  ) 

15.5 

21.7 

25.1 

Initial 
Delamination 
Stresses 
(  ksi  ) 

Gages 

20.4 

23.7 

22.1 

Silver  Ink 

17.2 

24.6 

24  9 

Acoustic 

Emission 

19.2 

23.2 

26.4 

Visual 

18.1 

23.5 

25.9 

Average  Stress  Calculated 
from  Continuous  traction  Q-23 
element  (ksi) 

15.883 

25.817 

29.935 
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632,  Delamination  of  [(0/-0)2/9O]3  Laminate 

A  sequence  of  tests  had  been  conducted  [5l]  to  monitor  the  material  damage  in 
[(0/— O)j/90],  laminate  specimens  under  incremental  loading.  The  value  of  0  were 

5°,  15°,  25°,  35°,  45°.  The  material  considered  here  was  T300/5208  graphite  epoxy 
with  the  following  elastic  constants; 

i:,,=22x]()'’i-)si 

I:,,=1.54xU)*’psi 

The  thickne.ss  of  the  ten-ply  laminate  averaged  around  0.06  inch  with  width  equal  to 
1  inch.  All  these  laminate  specimens  had  also  been  analyzed  [51]  using  the  assumed 
stress  hybrid  finite  element  model  [52]  in  conjunction  with  the  quadratic  tensor 
polynomial  failure  criterion  [53]  to  predict  the  onset  of  transverse  cracking  and 
delamination.  In  the  present  study,  continuous  traction  Q-23  element  with  a  uniform 
100  element-model  shown  in  Figure  (60)  was  used  to  analyze  one  quadrant  of  a 
typical  x=constant  plane  in  the  laminate  specimens.  The  calculated  stresses  along  the 
traction-free  edge  were  then  substituted  into  the  following  failure  criteria  to  determine 
the  possible  sites  for  initiation  of  transvers*-  cracks  and  delumination. 
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Figure  60:  KXFElement  Model 


6.3J2A  Anisotropic  Strength  and  Failure  Criteria 


With  macroscopically  homogeneous  but  orthotropic  materials,  development  of  a 
strength  theory  has  been  frequently  accomplished  by  extending  one  of  the  isotropic 
analyses  to  account  for  anisotropy.  Since  strength  theories  are  used  primarily  to 
predict  onset,  rather  than  mode  of  failure,  the  macroscopic  viewpoint  will  predominate. 
It  has  been  stated  [53]  that  all  the  failure  criteria  are  the  degenerate  cases  of  the 
tensor  polynomial  failure  criterion 


r  cr  +1'.  cr  cr  +F  .  cr  cr  cr,  +....^  1 

II  IJ  I  J  ijk  I  j  k 


(194) 


or,  explicitly 

Here,  cr,  are  the  stress  tensor  components  in  the  material  coordinates  and  F,,  F,j  and 
etc.  are  the  components  of  strength  tensors,  all  components  are  referred  to  the  material 
principle  axes.  In  (195),  terms  associated  with  cr^  ct^,  and  cr^  which  are 

F^,  Fj,  and  F^,  are  taken  to  be  zero  since  shear  strengths  are  the  same  for  positive  and 
negative  shear  stress.  It  is  also  assumed  that  there  is  no  interaction  between  shear 
stresses  and  normal  stresses,  thus  F^,  F^^,  F,^  etc.  become  zero. 

The  strength  and  failure  criteria  considered  in  the  present  study  include  the 
maximum  stress  criterion,  maximum  strain  criterion,  Hoffman's  criterion  and  the 
Tsai-Wu  criterion.  The  reason  for  choosing  these  criteria  was  not  only  their  popularity 
but  also  because  they  include  unequal  tensile  and  compressive  failure  strengths. 
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Maximum  stress  criterion 


I'ailure  of  material  is  assumed  to  occur  if  any  one  of  the  following  conditions  is 
satisfied  [47] 

o-  >X^  cr.>Y-;  cr  >7^ 

(196) 

cr^>R;  o‘j>S;  o'(,>T 

where  cr,,  cr,,  Cj  are  the  normal  stress  components;  O’^,  Cy  cr^  are  the  shear  stress 
comix)nenis:  \,,  Z,  are  the  lamina  normal  strengths  in  the  x,  y,  z  directions 
respectively;  and  R,  S,  'f  are  the  shear  strengths  in  the  yz,  xz  and  xy  planes, 
res|-)ecti\ely.  When  cr,,  ,cr^,  cr^  are  compressive,  they  should  be  compared  with 
'll  ,  and  Z, ,  the  normal  strengths  in  compression  in  the  x,  y  and  z  directions. 
Reddy  [54]  stated  that  the  maximum  .stress  criterion  could  also  be  expressed  in  the 
form  of  tensor  polynomial  criterion  as 

(cTj-X^Xo-i+X^Xo-j-YTXcr^+Y^Xa^-Z^Xo-j+ZpXcr^-RXo-^+R) 
(crj-SXcr^+SXo-j_-TX<r^+T)=0  (197) 

Comparing  (197)  with  (195)  and  ignoring  those  higher  order  terms,  the  strength  tensors 
are  [54] 


P  ^  J__J_  ,•  =_L__L  r  =  J _ -L 

■  ^  >^c’  ^  Y.^  Yp’  ^  Z,  Z,. 

p  = _ ] _ •  l'  = _ 1 _ •  F  = — ^ _ 


T.'r 


T  r 


F  —  ^  .  F  _  ^  «  F  _  ^ 

44  ^2’  5S  ^2-  .C  ^2 

P  _  I^.Fz.p  _  F.F3.P  _  ^3 

- 2"’  *^'3 - 2~’  - 2“ 


and  the  remaining  strength  constants  are  zero. 


(198) 
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Maximum  strain  criterion 


Failure  is  assumed  to  occur  if  one  of  the  following  conditions  is  satisfied  [48] 


€j>X^tI 

e,>R^;  €5>S,:€,>T^ 


(199) 


where  €,,  €3  are  the  normal  tensile  strains  in  the  x,  y,  z  directions  respectively: 

6^,  €5,  are  the  shear  strains  in  the  yz,  xz  and  xy  planes  respectively;  X,,,  Y,,,  Z,, 
are  the  tensile  strain  strengths  in  the  x,  y.  z  directions  and  R,.  S,.  T,  are  the  shear 


strain  strengths  in  the  yz,  xz  and  xy  planes  resjTectively.  Again,  expressing  the 
criterion  in  the  form  of  tensor  ixtiynomial  (195), 

(6  -X,,X€,+\^,Xe-Y,,,X€_,+Y^,.Xe-Z,,.X€3+Z,,Xe-R^Xe,+R,) 


(€ .-S  X€,+S  Xe, -T  Xe.+T  )=0 

€0  i 


(2(K)) 


Expressing  strains  in  terms  of  stresses  via  the  compliance  matrix  for  orthotropic 

materials,  (2(X))  can  be  expressed  in  the  form  of  (195)  and  we  have  [54] 

F  =F^+£i-LF''+-^F‘' 

1  c  2  s;  3 

^22  ^33 

F  -  ^‘,2  ^23  pA 

'2  c  12c  3 

^11  ^33 

F  =^F>^f!+fC 
^  S,,  '  S,,  ^ 


r'  1  j/^12'^2  1  i/'^13'i2  1  ^13  ^12  ^12  13 

^22  't’^C  ^33  ^33  ^22  ^22^33 

S. 


F  = 

Y.Y 


,  ^  *^23  ^2  1  ^12  y'Ai'^A  ^23  ^A^^A  ^12^23 

^33  ^11  *33  *11*33 


S.,S, 


T  C 


S ,  I 


F,  = -i- +( -^  )*—!—+( 

z  z  S  \  X  s 

T  <  '^11  't  r  ^22  'j'< 


8  8  8  8 

^12  rrA^^A  '^23  cApA  '^13  ,;j  p". 


'I  }^22 
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_  ^12  1  ^  ^^13^23  1 

‘2  5^3  ZtZc 


O^Q  Q  t'-c  c  ^c  iVc  c  ^c  -^2^3 


S.,S„  S. 


2  S,,S22 


^  ^1 1^33  ^33 


-  _  ^13  1  .  ^13  1  ,  ^12^23  _ 1 


S.,  S„  Z,..Z, 


11  ^33  S,,  ‘T'r 


.2  Y  \ 

V  •  -I-  •  /■ 


_l(_li^  +  l)i:V_l(Mli  +  4^)i:Vt-l(il4^  +  ^)rV^ 

TCC  >3')CC  C  I2')CC  c  23 

Z.  0,,033  ^  '^ll'^22  22  22  33  22 


p  _  ^23  1  -(-^^3  1  p  ^12^13  _ ^ 


S22  Y^-Yj,  S33  Z^Z(,  ^2^ 


_!(■  ^23  +i)pApA_lr^l2^23.  ^13  3pApA_  1/^23^13  .^12ApApA 

1  <1  <;  ^23  2  <;  s  s  12  O  S  s  s  13 

^  ^22^33  ^  ^11*22  ^11  ^  ^U'^33  ^11 


S22S33 

^33 

^13'*^  2  3 

1- 

c/l 

to 

^22^33 

^23^13 

(201) 


Here,  S,j  (i,j=l,2,3)  are  the  components  of  the  compliance  matrix,  and  Ff,  Fj,  Ff  are 
the  expressions  given  for  F,  F^,  Fj  in  the  maximum  stress  criterion. 


Hoffman's  criterion 

Hoffman's  criterion  [49]  is  a  special  case  of  the  tensor  polynomial  criterion  for  the 
following  choice  of  the  parameters  F,  and  F,j: 
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..  1  1 


1  1  ,,  1  i 


•  '  x^.  z,,  z,. 

F  =_I— F  =-1— F  =-!- 

“  X^X^’  22  Y^Yp’  Z^Zj, 

F..=-p:  ’’“'P 


■  ( 

2 

•S  y  - 

T  C 

Vc 

’(  ’ 

2  \.,X, 

.+  ’  _ 

1 

z,z. 

Y,Y, 

2  7  7 

+  ’  _ 

1 

Y.,x;. 

X.,X,. 

Tsai-Wu  criteria 

The  Tsai-\Vu  criterion  is  given  by 


where 


F  cr  +F  cr  cr  ^  1 

II  ij  I  j 


p  _ L.  F  =J _ L.  F  =j _ L 

I  XF  >  Y,  YF  n  Zc 

tr  —  ^  ,1::^  ^  ,1:;^  ^ 

^‘1  XX  '^22  yY  ’^33”7  7 
A^C 

P  —  ^  •  p  —  ^  .  p  _  ^ 

^4  s.  j.a’^e  ^a 


»^a=- 


^.t=- 


2Vx,x^y,v,. 


2  Y 


Fa3=-: 


2VY^Y^Z,Z^ 

Here,  it  is  noted  that  the  maximum  stress  and  maximum  strain  criteria  involved 
several  separate  equations,  and  there  was  no  allowance  for  interaction  of  the  stresses  or 
failure  modes.  The  Hoffman  criterion  and  Tsai-Wu  criterion,  however,  do  provide  for 
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interaction,  and  the  interaction  is  fixed.  That  is,  these  failure  expressions  are  not 
invariant  with  respect  to  coordinate  system.  .As  expressed  in  terms  of  quadratic  tensor 
polynomial  shown  in  (202)  and  (204)  respectively,  the  only  difference  between  these 
two  criteria  was  on  the  determination  of  strength  tensors  Fj^,  F,j,  F^j. 

Onset  of  Material  Damage 

I'.ach  laminate  specimen  was  tested  individually  in  an  electro-hydraulic, 
ser\o-con trolled  closed-loop  testing  machine  [5l].  fhe  strain  and  nominal  stress  at  the 
first  sight  of  transverse  cracking  and  onset  of  delamination  are  summarized  in  'I'able 
(5). 

The  measured  strength  (ksi)  of  'O(K)/5208  graphite  epoxy  are  given  by  [50] 

Longitudinal  tension  :  X.j.  =  210 
Longitudinal  compression  :  =  200 

Transverse  tension  :  Y.j,  =  10 

Transverse  compression  :  =  21 

Shear  in  1-2  plane  :  S  =  13 

It  is  further  assumed  that  Z(.=Yc,  R=S  and  T=S/2.  Substituting  above 

information  into  (198),  (201),  (202)  and  (209)  to  calculate  Fj  and  F,j,  the  strength 
tensor  for  any  complex  stress  state  can  then  be  obtained  and  compared  with  the  actual 
stress  tensor.  Failure  is  a.ssumed  to  occur  when  the  magnitude  of  the  actual  stress 
tensor  exceeds  that  of  the  strenth  tensor. 

Transverse  Cracking 

Based  on  the  stress  field  calculated  from  the  continuous  traction  Q-23  element,  the 
four  failure  criteria  discussed  in  6.3.2.1  were  applied  to  every  point  along  the 
traction-free  edge  for  all  five  laminate  specimens  at  the  strain  levels  shown  in  Table 
(5),  respectively,  when  the  first  sight  of  transverse  cracking  was  detected.  The  results 
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I.lo  1.00  2.00  3.00  11.00  5.0 
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15) 


gure  63:  Ratio  of  Slress-to- Strength  Tensor  for  [(±25)2 

under  Transverse  Cracks  Loading 
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Figure  64: 


Ratio  of  Stress-to-Strength  Tensor  for  [(±35)^  Laminate 

under  Transverse  Cracks  Loading 


Z/H 

t.M  1.00  2. 70  3.00  1..  00  5.00  6.00  7.00  8^.00 _ 9^.0^ 


-1.00  0.00  I.OQ  2.00  3.00  ^.00  5.00  6.00  7.00 

Figure  65:  Ratio  of  Stre.ss-to-Strength  Tensor  for  [(±45)2 

under  Transverse  Cracks  l.oading 
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ure  shown  in  Figures  ((>1  )-(65)  along  with  Chou's  [30]  predictions.  Those  points  tliat 
lie  in  the  region  where  the  actual  stres.s-to-strengtli  tensor  ratio  is  greater  than  unity 
represent  failure.  Due  to  the  discrepancy  of  the  calculated  stress  field  based  on 
different  numerical  schemes,  the  prediction  of  the  lamina  failure  surface  from  the  same 
failure  criterion  (such  as  Tsau-Wu  theory)  varied  significantly  through  the  laminate 
free-edge.  An  obvious  failure  phenomenon  resulting  from  the  transverse  cracking 

within  the  9()-degree  layer  was  delected  based  on  the  continuous  traction  Q-2.1  element 
li>r  all  the  laminate  specimens.  At  the  same  time,  initiation  of  edge  delamination  at 

the  interlaces  lietween  H  and  —0  become  apparent  as  the  \alue  of  0  increases,  winch 
was  not  indicated  according  to  ('hou's  analysis.  However,  the  fact  that  transxerse 
cracks  always  tKcurred  prior  to  delaminaiion  in  all  cases  is  noticed,  and  this  indeed 

matches  experimental  observation  [5l].  Here,  it  is  noted  that  the  magnitudes  of 
stress- to-strength  ratios  shown  in  these  figures  sometimes  departed  significantly  from 
unity  particularly  in  the  interior  of  transverse  layer  and  near  the  interfaces.  This 
could  possibly  be  due  to  the  inaccurate  insitu  transverse  strength  data  and  the 

inappropriate  assumption  of  the  interlaminar  strengths. 

Edge  Delamination 

Following  ihe  same  pr(x;edure  as  in  the  prediction  of  transverse  cracking,  maximum 
stress,  maximum  strain,  1  lol  l  man  and  Fsai-Wu  criteria  were  applied  to  every  point 
along  the  I'ree-edge  of  various  laminate  specimens  at  the  respective  strain  level 
correponding  to  the  onset  of  delamination.  For  illustration,  failure  surfaces  predicted 
from  the  continuous  traction  Q-23  element  for  0=5°  and  0=25°  laminate  specimens  are 
shown  in  Figures  (66)  and  (67).  In  the  case  of  0=5°,  Figure  (66)  shows  that 
following  the  transverse  cracks  formed  in  the  90-degree  layer,  delaminations  were 
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ileveloi')ed  al  the  interfaces  between  5/-5  layers.  Meanwhile,  transverse  cracks  also 
extended  to  the  angle-ply  layers  under  incremental  loading.  The  fact  that  all  the 
failure  surfaces  exceeded  unity  shown  in  Figure  (67)  might  result  from  the  inaccurate 
material  strength  data.  For  the  6=25°  laminate,  however,  transverse  cracks  were  still 
confined  to  the  90-degree  layer  as  delamination  propagated  at  the  25/-25  interfaces. 
Rased  on  these  observations,  we  can  infer  that  fiber  orientations  of  the  laminates  with 
tile  same  stacking  set|uence  iia\e  played  an  impintant  role  on  the  determination  of 
damage  modes  under  incremental  hxiding. 

In  general,  the  llolfman  theory  had  more  conser\ati\e  predictions  than  the  others 
on  the  initiation  of  transverse  cracking  within  the  9()-degree  layer,  and  the  maximum 
strain  criterion  predicted  conservatively  on  the  subsequent  edge  delamination  at  the 
interfaces  between  angle-ply  laminae.  Since  the  materials  were  assumed  linear  elastic 

in  the  analysis,  the  applied  strain  loading  corresponds  to  the  onset  of  transverse  cracks 
or  delamination  based  on  the  continuous  traction  0-23  model,  and  the  failure  criteria 
would  be  expected  to  be  lower  than  the  experimental  observation.  However, 
throughout  the  analysi.s,  delamination  was  assumed  to  occur  in  a  state  of  generalized 
plane  strain  without  the  influence  of  tran.sverse  cracking.  In  reality,  this  is  not  the 

case.  More  work  needs  to  lie  done  to  study  the  interrelationships  between  delamination 

and  other  damage  mcxles  such  as  matrix  cracking  and  fiber  breakage,  etc.  Also,  many 
practical  composite  systems  actually  exhibit  extensive  nonlinear  mechanical  response  in 
shear  and  transverse  to  the  reinforcement,  re.sulting  in  nonlinear  laminate  mechanical 
behavior.  Extension  of  the  present  continous  traction  finite  element  procedure  to 
include  nonlinear  material  behavior,  along  with  careful  determination  of  material 

properties  and  strength  data,  may  lead  to  better  estimation  of  initiation  of  various 
damage  modes  under  incremental  loading. 

l.-Sb 
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Figure  66: 


Ratio  of  Stress-to- Strength  Tensor  for  [(±5)2  Laminate 

under  Delamination  Loading 
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SECTION  VII 


DISCUSSION 

The  problem  of  I'ree-edge  cielamination  in  composite  laminate  coupons  subjected  to 
unilorm  inplane  estension  has  Ireen  incest igateii.  llelore  delaminaiion  can  be  predicted 
\  i.i  ;i  stress  based  laiUire  criterion,  an  acciir.ite  stress  c.ilculation  within  the  laminate, 
particuhirlv  netir  the  interlaces  and  tiiiction  free  boundary,  is  necessary.  However,  the 
stress  field  under  this  sitiiaiutn  is  highly  comple.x  in  nature.  Besides,  the  anisotropy 
and  heterogeneity  of  the  material  system,  and  presence  of  the  traction-free  boundary 
makes  the  analysis  difficult.  Literature  on  this  subject  w’as  abundant  but  an  effective 
as  well  as  reliable  solution  had  not  been  found. 

The  present  research  effort  has  resulted  in  developement  of  a  continuous  strain 
finite  element  model  in  plane  elasticity  based  on  the  compatible  cubic  interpolation 
function  proposed  by  ('lough  and  I'elippa  [40],  in  which  the  normal  slope  continuity 
was  ensured  across  the  interelement  boundary.  Lxtending  the  continuous  strain  model 
to  tlie  analysis  ol  a  pseudo  two-dimensional  free-edge  delamination  coupon  under 
unilorm  extension,  a  continuous  strum  field  along  both  inplane  and  transverse  directions 
was  obtained.  However,  due  to  material  anisotropy,  the  stresses  along  the  interfaces 
between  differently  oriented  layers  w'ere  discontinuous.  Also,  traction-free  boundary 
was  not  satisfied.  The  continuous  strain  model  w'as  used  as  the  basis  for  the 
developement  of  a  continuous  traction  finite  element  procedure.  Knowing  the  fact  that 
the  displacement  field  within  each  element  is  destiihed  by  nodal  displacement 


components  ami  their  gradients,  t('  ensure  traction  continuity,  a  transi'ormation  procedure 
was  developed  to  map  the  gradients  normal  t(r  element  Ixtundary  to  a  mixed  set  ol 
degrees  of  freedom  through  appropriate  displacement-stress  relationships.  For  global 
assembly,  the  nodal  degrees  of  freedom  of  this  element  include  interlaminar  stress 

components  at  the  corner  nodes,  as  well  as  traction  components  at  the  mid-side  nodes 
of  each  element.  This  ensures  continuity  oi'  displacement  and  traction  along 
interelement  boundaries  as  well  as  across  laminate  interfaces  providing  a  small 
delormaiion  situation  is  considered.  At  the  same  time,  equilibrium  condition  is 
maintained  between  twii  adjacent  elements  (layers).  A  significant  aspect  of  this 

displacement-based  formulation  is  that  it  allows  tracticm-free  boundary  conditions  to  be 
sfiecified  in  a  point-wi.se  sen.se. 

In  the  four-ply  laminate  analysis,  numerical  results  calculated  from  the  continuous 
traction  Q-23  element  generally  agreed  well  with  Pagano's  analytical  solutions  [8] 
although  these  two  schemes  were  based  upon  quite  different  theories.  For  illustration. 
Table  (6)  outlines  the  basic  characteristics  associated  with  each  of  these  approaches. 

The  approximate  solutions  for  stress  components  and  cr^,  which  play  an  important 
role  in  delamination  of  composite  laminates,  were  calculated  using  both  approaches  and 
found  to  have  similar  distribution.  'I'he  study  also  revealed  that  the  pattern  of  mesh 
refinement  liad  significant  effect  on  the  estimates  of  interlaminar  stress  field  in  the 
vicinity  of  traction-free  edge  or  near  the  interface  betw'een  two  differently  oriented 
layers.  Here,  it  is  essential  to  realize  that  the  continuous  traction  finite  element 
procedure  is  only  applicable  to  the  0-23  element  and  cannot  be  simplified  to  0'f5  and 
0-19  elements.  This  is  because  the  continuity  of  traction  across  laminate  interface 
cannot  be  simplified  in  the  absence  of  mid-side  nodes  at  the  interface  between  two 
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Table  6;  ('omparison  ol  Pagano's  ’I'lieory  and  (^oniinuous  'I’raction  Q-23  Hiement 


Method 

Pagano 

Continuous  traction  Q-23 

Variational 

principle 

Reissner’s  variational 
principle 

Minimum  potential 
energy  principle 

Type  of  formulation 

Mixed 

Displacement 

Basis  of  field 
equations 

Plate  theory 

Elastic  solid 

Assumed  inside  each 
layer  (element) 

Kinematic  relations 

Ic  Stress  equilibrium 

Kinematic  relations 
&  constitutive  lav; 

Along  interlaminar 
boundary 

Continuous  traction 
t  weighted  displacement 

Continuous  displacement 
&  traction 

Assumed  stress 
v'.r.t.  Z-a::is 

Inplane - linear 

Transverse - cubic 

Quadratic 

Unknovms  in  final 
equations 

Weighted  displacements 
&  interlaminar  stresses 

Displacements 
&  interlaminar  stresses 

Solution  technique 

Direct  solving 
differential  equations 

Finite  element  method 
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adjilcent  layers.  In  cdtnparison  with  the  continuous  strain  y-23  element,  the 
intrcxluction  ol'  the  translormation  prttcess  in  the  element  makes  the  continuous 

traction  procedure  more  expensive.  However,  the  continuous  traction  Q-23  element 
significantly  improves  the  reliability  of  the  stress  field  solution  because  of  the 
interlaminar  stress  continuity  at  the  interface  between  differently  oriented  fiber  layers, 
along  with  satisfaction  of  traction-free  boundary  condition  along  specimen  edges. 

Ajiplication  to  the  multiply  laminate  six'cimens  with  stacking  sequence  of 
[id 01^.,  bO  J  further  illustrated  the  potential  (>1  the  ei'ntinuous  traction  finite  element 
procetiure  in  the  analysis  (if  edge  delaminalion  problem.  Satisfactory  agreement  was 
generally  ohser\ed  for  interlaminar  sire.ss  distributions  as  well  as  laminate  displacement 
field  between  continuous  traction  f,)-23  element  and  constant  strain  element  solutions 
[44]  except  that,  in  the  vicinity  of  the  free-edge,  the  constant  strain  element  w’as 
deficient  due  to  nonsalisf action  of  traction-free  boundary  condition  (Ty,=0)  and  the 
assumption  of  (r„=0)  imbedded  in  the  axisymmetric  analysis.  The  results  from  the 
continuous  traction  Q-23  element  would  be  expected  to  be  sujjerior  to  the  constant 
strain  element  (conventional  assumed  displacement  elements)  for  prediction  of  stress 
field  in  the  free-edge  delamination  specimens.  Of  course  the  simple  axisymmetric 
mixlel  is  economical  to  use. 

Regarding  prediction  ol  damage  initiation  in  laminate  composite  couperns, 
under  incremental  loading,  the  continuous  traction  finite  element  procedure 
along  with  some  popular  anisotropic  failure  criteria  was  found  to  be  successful  in 
modelling  some  failure  phenomena  observed  in  the  experiments.  Basically,  the  laminate 
specimens  under  analysis  were  assumed  to  be  made  of  linear  elastic  brittle  materials. 
Thus,  initiation  of  delainination  directly  led  to  catastrophic  laminate  failure  regardless 
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of  the  damage  accumulation  prcKess.  Numerical  e\i-)eriment  discussed  in  the  last  section 
revealed  that  the  lloffman  criterion  had  a  more  conservative  prediction  on  the 


initiation  of  transverse  cracking  within  the  90-degree  layer,  and  the  maximum  strain 
criterion  on  the  subsequent  edge  delamination  between  the  angle-ply  laminae  interfaces. 

In  summary,  we  conclude  that  the  proposed  continuous  traction  finite  element 
scheme  not  only  overcomes  the  drawback  of  deficient  stress  calculation  arising  in  the 
conventional  assumed  displacemeni  approacli.  hut  also  provides  a  reliable  as  well  as 
erie>.ti\e  numerical  solution  procedure  with  a  wider  range  of  applicability  to  the 
analvsis  of  the  I  ree-edge  delamin.ition  jirohlem.  fhough  based  on  a  completels 
different  variational  I'ormulation,  this  model  has  shared  the  characteristic  of  continuous 
displacement  as  well  as  traction  fields  across  laminate  interface,  with  Pagano's 
approximate  theory  derived  from  a  mixed  formulation.  Though  developed  for 
evaluation  of  stresses  in  composite  laminates,  the  continuous  traction  procedure  is  also 
applicable  to  analysis  of  layered  media  involving  material  interfaces  where  a 
two-dimensional  or  pseudo  two-dimensional  representation  is  applicable.  This  would 
include  stresses  in  layered  airfield  and  highway  pavements,  pressures  on  tunnel  lining, 
etc. 
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Appendix  A 

DERIVATION  OF  COMPATIBLE  CUBIC 
INTERPOLATION  FUNCTIONS 


This  <ippendi\  conuiiiis  a  summary  ('I  I  eli]-)pa's  141.42]  approach  in  deriving  tlie 
cubic  compatilde  interpolations  for  in-plane  displacement  u  in  a  more  detailed  i'ormat. 

In  order  to  derive  the  cubic  interjxdation  1  unctions  for  the  complete  triancular 
elentent,  three  different  coordinate  systems,  i.e.  triangular  coordinate,  local  and  global 
Cartesian  ctxtrdinaies  should  be  defined  as  illustrated  in  Figure  (A.l).  The  geometry  of 
an  arbitrary  triangular  element  can  be  expressed  in  a  Cartesian  coordinate  system  by 
its  nodal  coordinates  or  its  projected  dimensions  as  shown  in  Figure  (A.2),  or 
alternately  by  its  intrinsic  dimensions  as  defined  in  Figure  (A.3). 

Let  i  and  k  denote  the  first  and  .second  cyclic  permutations  of  i=l,2,3  (i.e.  j=2,3,] 
and  k  =  3.1,2),  the  projected  dimensions  may  be  written  as 

a  =  X  -  X  :  b  =  V  -  \’  (A.l ) 

'  k  I  I  j  ■  k 

Also,  the  intrinsic  dimensions  may  be  defined  in  terms  of  the  projected  dimensions. 
Referring  to  Figure  (A.2),  define 

d 

A  =  —  (A.2) 

'  1 
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X 


(Kill  and  (iiobal  ('artesian  ('oordinates 


/jt  —  ]  —  K 

'  1  \ 


(AJ) 


(a  u.  +b  b, ) 

d  =  lJ; - (A.4) 

‘  ] 

1 

'I'he  triangular  coordinates  4,,  of  any  point  ”P”  in  the  triangle  inay  be  defined 

eitlter  as  the  ratios  of  the  areas  A,  ol'  the  suhinangles  subtended  by  the  point  to  the 
total  area  A  of  the  triangle,  or  as  the  ratios  t'f  the  norma!  distances  n,  to  the  Iteight 
h,.  i.e. 


(A. 5) 


as  shov^n  in  Figure  (A.l).  It  is  noted  that  the  triangular  coordinates  are  related  by 
the  constraining  condition 

With  reference  to  Figure  (2).  the  displacement  interpolation  functions  for  each 
subelement  (i)  express  the  relationship  between  the  displacement  u'"  within  the  element 
and  the  ten  displacement  components  of  its  nodal  jxiints  r''  as  follows 


u'"  =  {0‘'V{i 


I'or  example,  the  nodal  displacement  vector  for  subelement  1  is 


(A. 6) 


<r“V={u^,u, 


,U,  ,u^,u^  ,u^,  ,u^  ,u_ 


'3  >3 


(A.l) 


The  corresponding  ten  cubic  interpolation  lunctions  expressed  in  triangular  coordiates  are 
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</a,  4,-a,  i,)+(a^  /i,  -a, 

^3~^3  "‘^3  ^3 

^J(3-2{p+6Xy'<,i,$3 

<^(3-2^3) 

^3<^'''‘^.-^>V’^2> 

i^^b\%-b:\^) 


where  the  suhsenpts  correspKtnd  to  the  renumbered  nodes  ol  the  subelemeni.  and  are 
the  local  triangular  coordinates  of  points  in  subtriangle  1.  With  this  convention,  the 
interpxDlation  functions  for  subelements  2  and  3  are  the  same  as  (A.8)  appropriately 
permuting  the  subscripts  and  superscripts.  It  should  be  noted,  how'ever,  that  the  nodal 
displacements  in  (A.b)  are  identified  by  node  numbers  defined  for  the  complete  element 
assembly. 

If  the  vector  (f)  of  all  ncxlal  displacements  of  the  complete  element  assembly  is 
written  as  the  ordered  set 


{r}*={u.,u  ,u  ,u,.u  ,u  ,u,,u  ,u  .u  ,u  ,u  lu  ,u  ,u  ) 
'  ''1  >1  ’  '2  >2  3  '3  ''3  "d  'S  %  V  >0 


the  displacement  in  subelement  1  can  be  expressed  as 


(p  j  I  1 1  - 1  r^‘^'  1 

u  =  1$  )  Irl  =\<t)  \<f>  J 


(A.IO) 
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where  {$“'}  is  similar  to  (A.8J,  bui  expanded  witli  5  zeros  to  account  lor  the  nodal 
displacements  not  assixjiated  with  subelement  1,  and  with  appropriate  arrangement  ol 
terms.  The  vectors  {c{>l"}  and  represent  the  interpolation  functions  for  the 

external  and  internal  nodal  displacements  respectively. 

Hxpressing  the  displacements  in  the  other  subelements  similarly,  the  complete 
system  oi  displacements  can  Ire  written  as 


•.  1 

0" 

1  0;,' 

u 

V 

r 

.  ■»  1 
u  '  ' 

= 

•  0';"' 

— 

■  n 

r 

u 

1 

(.A.m 


(A.ll)  is  an  expression  of  the  cubic  displacement  patterns  in  the  three  subelements. 
The  displacement  of  two  adjacent  subtriangles  are  identical  along  their  common 
boundary.  The  normal  slofie  at  any  of  these  nodes  (say  7  of  subelement  1)  is  given 
by 


r 


( )  (  )  90  90 

where  b",.,  b',',  respectivelv  are  values  of  { -  -  at  node  7  for  subelement  1, 

9”  9n 

and  n  denotes  the  axis  normal  to  the  element  boundary.  To  maintain  internal  slope 
continuity,  it  is  necessary  that  where  the  negative  sign  results  from  the 

convention  that  the  positixe  normal  is  directed  outwards.  I’or  the  three  points  7,  8,  9, 
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(A. 13) 


bV>b;^’ 

,  (1)  .  ,  (3)1 
^7o+No 

r 

0 

e'; 

= 

bf+b''^ 

— 

0 

e(^3)+0(^2) 

,  (3)  ,  .  (2) 
ba  +b. 
9d  9o 

r 

o 

0 

or  symbolically 


Ir 


=  0 


(A. 14! 


The  values  of  r„  which  will  satisfy  these  compaiibility  conditions  are  obtained  by 
solving  (A. 14),  i.e. 


r  =-B  'Br  =  Lr 

^  O 


(A.15) 


Substituting  the  slope  continuity  constraint  of  (A.15)  into  (A.ll),  the  fully  compatible 
displacement  field  in  the  three  subelements  becomes 


• 

(1) 

u 

(2) 

U 

= 

+ 

,(2'' 

0 
^  1' 

L  |r)  = 

* 

(.1) 

u 

ir} 


Explicit  expressions  of  these 


functions  for  subelement  3  are 


(A. 16) 
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=  r>3^,-a,<3)+(a3M3-a.)^  ,4,^3 


+lCj3(2a  -a3/X3-a2X3)<,+3(a -a3M3X2+(-3a,+2a3M3+a3X^){3] 


>  I 


+  lU'l3(-2bj+b^M,+b^X;)^^  +  3(b^M  -bj)^.+(3b  -2b,M  -b 


4>;^''  =  r;(3-2^3)+6X,^,^,^,+r}3(M'-X3)^3+(2X,-/.,)^3-,U3<,] 


^  1^3-S^  I  ^+^^2-“3^3^^  1^3 

\  -> 


+  l5j3(a3X3-a3H,+3(a,M,+a3X3-2a2)<3+(3a3-a,M,-2a3^3^<3J 


^u''’  =  <>3^,-^^3)+(b3X3-^2)<,<2^3 

y  0 


+  lr’[3(b,-b,X,H,  +  3(-b,M,-b3X3+2b,H,+(-3b3+b,M,  +  2b3X3)^3] 


<t.;''  =  4;i3(l+/u,H,+3(1+X,)^,+(l-M,-X,)iJ 


^ =  Tri(2( a ,  +  3a  _,+a  , -3(  3a ,  +a  3+a  ,  X ,  )^  ,+(a ,  X ,  -a  3] 

x,  o 


J76 


|■^5[-3(b,+3b^+b,/xp^,+3(3b,+b^+b,\^)^,-(b,X 

^3 


"4  -^^3 

"5  '’*  1 

”(>  ‘^2 

The  above  set  of  interpolants  is  applicable  to  all  points  lying  in  subtriangle  3.  For 
points  lying  in  subtriangle  1  and  2,  {4>‘’^},  can  be  written  down  by  cyclic 

permutation  of  all  subscripts  and  superscripts  in  (A.17).  All  the  symbols  on  the  right 
side  of  (A.17)  relate  to  the  complete  triangle. 
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LIST  OF  SYMBOLS 


A  list  of  the  most  commonly  used  symbols  and  their  general  meaning  follows. 


a,  b 

A 

A 

1 

M 

C, 


s 

IJ 


d 

1 

'^1 


i: 


Global  dimensions  of  a  triangle 
Area  of  triangle 
Area  of  subtriangle  i 
Material  constants 

Components  of  the  stiffness  matrix  in  the 
global  coordinate  system 

Components  of  the  stiffness  matrix  in  the 
material  coordinate  system 

Components  of  the  compliance  matrix  in  the 
global  coordinate  system 

Components  of  the  compliance  matrix  in  the 
material  coordinate  system 

('omponents  of  the  translormed  reduced 
stiffness  matrix 

(Components  of  the  reduced  still  ness  matrix 
Projection  of  a  corner  of  a  triangle  over  opposite  side 
Inversion  of  O 

•j 

('omponents  of  the  reduced  material 
compliance  matrix 

'foung's  moduli 
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1 


('artesian  eoniimnenls  nl  tlie  lx>dv  lorce  vector 


I. 


h 

K 

IJ 

1 

1 

\ 

A 

(' 

I) 

R 

S 

Sp 

R 

n, 

QR 

a 

\\ 

u 

1 


Triangle  heights 

lamina  thickness 
I'lements  ol  Stiffness  matrix 

'I'riangle  side  lengths 
Normal  tci  boundary 

Linear  oirera'or  or  matrix  ol  linear  o]X'rators 
on  a  region  R 

Linear  operator  or  matrix  oT  line.ir  oirerators 
on  the  IxHindarv  ol  R 

Domain  ol  operator  A 
n-dimensionai  Luciidean  space 

Open  connected  region  in  ii" 

Boundary  of  R 
(Complementary  subset  of  S 

Closure  of  R 

L(Ku1  cartesian  components  of  the  unit 
normal  to  a  surface 

Local  cartesian  components  of  the  unit 
tangential  to  a  surface 

Open  connected  subregion  of  R 

Boundary  of  subregion 

('artesian  conifxments  of  the  stress  tensor 

(Cartesian  components  of  the  displacement  vector 

('artesian  components  I'l  tin  traction  sector 
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{'artesian  comixments  ol  tlie  prescribed 
traction  vector 

Kronecker's  delta 

Cartesian  components  of  the  isothermal 
elasticity  tensor 

Interpolation  functions 

Assumed  displacements  in  element  m 

Assumed  i>encrali/.cd  nodal  displacen'ents  at  the 
Iroundary  of  element  m 

Assumed  j^eneraii/ed  nixlal  displacements  internal 
to  element  m 

Vector  of  strains  in  element  m  deriving  from  u"‘ 
Stiffness  matrix 
Load  vectors 

Cartesian  coordinates  in 
Poisson's  ratio 
Shear  mcxluli 

Components  of  infinitesimal  or  linear  strain  tensor 

Natural  coordinate 

Rotation  angle  from  global  axis  x 
to  material  axis  1 

Transformation  matrix 
Applied  uniform  strain  loading 

Displacement  component  in  x-direction 
Displacement  component  in  y-direction 
Displacement  component  in  z-direction 


ISO 


Displacement  function  in  x-direction 
Displacement  function  in  y-direction 
Displacement  function  in  z-direction 
fiber  orientation 
Geometric  parameters 
Linear  functional 

IDisplacement-stress  transformation  matrix  for  corner  node 
Displacement-stress  transformation  matrix  for  mid-side  node 
Bilinear  mapping  on 

Coordinate  transformation  tensor 
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